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1.0  OBJECTIVE 


The  general  objective  of  this  contract  is  the  development  of 
a computer  program  for  calculation  of  beam  trajectories  in 
coupled-cavity  traveling  wave  tubes. 

Specific  objectives  set  out  in  the  statement  of  work  include 
the  following; 

i)  tlie  program  will  be  in  FORTRAN  IV  level  H. 

ii)  the  beam  will  be  represented  by  a disc  model  up  to  the 
beginning  of  the  saturation  region. 

iii)  the  beam  will  be  represented  by  a ring  model  of  at  least 
96  rings  per  wavelengtti  in  the  saturation  region. 

iv)  the  speed  of  the  program  shall  allow  the  96  ring  calcu- 
lation to  be  carried  out  in  5 minutes  of  CPU  time,  or 
less,  per  cavity,  on  an  appropriate  computer. 

v)  che  program  shall  include  a self-contained  routine  to 
generate  an  rf  vector  potential  matrix,  to  avoid  depen- 
dence on  the  Los  Alamos  program  LALA. 

Though  not  stated,  it  was  imderstood  that  the  program  would 
include  rf,  magnetic  and  space  charge  fields,  without  paraxial 
approximations,  and  that  the  magnetic  fields  should  include 
both  uniform  (solenoid)  and  nonuniform  (PPM)  cases.  It  was 
also  understood  that  the  interaction  between  the  beam  and  the 
rf  fields  would  be  computed  in  both  directions  • — that  is,  the 
fields  would  be  appropriately  modified  by  the  computed  beam 
trajectories,  not  merely  applied  from  external  sources. 

This  report  describes  the  analytical  background  to  the  devel- 
opment of  the  computer  program.  It  is  not  necessary  to  read 
the  report  in  order  to  use  the  program:  a separate  User's 
Manual  gives  all  the  instructions  necessary  for  setting  up  a 
case  and  interpreting  the  results.  But  familiarity  with  this 
report  is  necessary  for  anyone  intending  to  modify  the  program 

The  analysis  is  specific  to  coupled-cavity  TWT  circuits  at 
this  stage,  but  much  of  it  is  sufficiently  general  for  future 
application  to  other  0-type  tubes. 


2.0  DESCRIPTION  OF  THE  PHYSICAL  MODEL 

2 . 1 The  Tube 


The  tube  will  be  represented  as  a sequence  of  gaps  and  tun- 
nels, as  shown  in  Figure  1;  there  are  rf  voltages  across  the 
gaps,  determined  by  the  rf  power  flowing  in  each  cavity,  and 
the  rf  fields  due  to  any  one  gap  are  taken  to  extend  into  the 
tunnels  on  either  side  as  far  as  the  midplanes.  Beyon.;  these 
planes  the  fields  due  to  the  adjacent  gaps  take  over.  This 
assumption  that  the  fields  due  to  one  gap  become  negligible 
beyond  the  midplanes  is,  of  course,  not  exact;  but  for  typical 
tube  structures  the  fields  at  these  planes  are  25  to  30  dB 
below  the  gap  fields,  so  that  it  is  a reasonable  simplifying 
assumption.  A numerical  example  supporting  this  will  be  foimd 
at  the  end  of  Section  4.7. 


2.1.1  Specific  Model 

In  order  to  have  a consistent  set  of  test  cases  for  numerical 
trials  and  illustrations,  an  imaginary  (but  not  imrealistic) 
tube  design  was  constructed. 

Taking  a goal  of  50  kW  peak  output,  30%  bandwidth  centered  on 
10  GHz,  a preliminary  rule-of-thumb  TWT  program  gave  36  kV, 

1.2  qp  for  the  beam,  .203"  for  the  turjiel  diameter,  .297"  for 
the  cavity  period,  9 ohms  interaction  impedance,  and  1.02x10° 
m/s  phase  velocity  (1.48tt  per  cavity)  at  10  GHz.  The  expected 
electronic  efficiency  was  24.9%. 

After  adjusting  the  voltage  upward  to  38  kV  at  1.1  pP  to  allow 
for  relativistic  effects  not  included  in  the  simple  program, 
and  rounding  off  other  parameters  to  convenient  values,  the 
following  set  of  nominal  parameters  was  adopted: 

Tube  type:  ’Navtest' 

Frequency;  10  GHz 
Power  output;  50  kW  peak 

Tunnel  diameter  .2”;  cavity  period  .3"»  magnet  period  .6" 
gap  .1". 

Beam  38  kV  1.1  pP  (approx.  8 amps),  b/a  = .7. 

Cavities  = 30;  impedance  10  ohms;  loss  0.1  dB/ cavity; 
sever  at  cavities  12  and  13,  phase  velocity 
1.0x10°  m/s. 

The  tube  structure  is  shown  in  Figure  2(a). 


Figure  1 


Idealized  TWT  Structure 


Figure  2(a) 

Structure  of  Coupled-Cavity  TWT  'Navtest' , 
used  in  numerical  examples 

Positions  of  current  loops  used  in  magnetic  field 
representation 
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This  case  was  run  on  the  large  signal  program  with  the 

results  shown  in  Figure  2(b).  The  upper  block  shows  the  input 
data,  followed  by  various  derived  quantities,  including  the 
equivalent  Pierce  parameters.  The  lower  block  shows  the  power 
saturating  at  61.8  kW,  l.e.  about  1 dB  margin,  at  43  dB  gain. 
The  energy  balance  in  the  last  column  is  within  .3  dB,  which 
is  quite  satisfactory.  The  output  is  plotted  in  Figure  2(c) 
showing  a very  normal  type  of  Applegate  diagram  for  a high 
power  over-voltaged  tube.  The  electronic  efficiency  of  22.2% 
is  somewhat  less  than  the  24.9%  estimated  by  the  preliminary 
program,  but  not  unreasonable.  Overall,  this  seems  to  be  a 
self-consistent  design  for  program  test  purposes,  and  its 
parameters  will  be  used  for  the  test  cases  for  the  rf  field, 
magnetic  field,  etc. 


2.2  The  Beam 

The  beam  will  be  represented  by  a one-wavelength  segment, 
traveling  down  the  tube  at  the  dc  beam  velocity.  The  asstunp- 
tion  is  made  that  this  wavelength  is  preceded  and  followed  by 
identical  wavelengths:  this  assumption  allows  us  to  do  two 
things : 

i)  compute  space  charge  forces  by  a fast  Fourier  analysis 
method,  which  implies  that  the  segment  considered  is 
part  of  an  infinite  sequence  of  identical  segments; 

ii)  replace  any  element  of  the  beam  which  leaves  the  segment 
at  one  end,  by  a corresponding  element  entering  at  the 
other  end;  i.e.,  an  element  can  always  be  moved  up  or 
down  one  beam  wavelength  to  keep  it  in  our  working  range. 

Since  the  tube  is  intended  to  be  an  amplifier,  the  bunching 
in  general  increases  along  the  tube,  so  the  assumption  of 
identical  wavelengths  ahead  and  behind  cannot  be  strictly 


[13  J.R.M.  Vau^an,  'Calculation  of  Coupled- Cavity  TVfT 
Performance',  IEEE  Transactions  on  Electron  Devices, 

ED-22  #10,  October  1975,  pp.  880-890. 

* References  v;ill  appear  as  footnotes  on  the  pages  where 
they  first  occur,  and  will  also  be  collected  in  a com- 
plete list  at  the  end  of  this  report. 

I See  page  93.  I 
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correct.  However,  the  tunnel  walls  exert  a shielding  effect 
which  diminishes  the  effect  of  more  distant  charges  very 
rapidly,  so  that  errors  in  estimating  their  magnitudes  have 
very  little  effect  on  the  final  results.  Indeed,  the  real 
reason  for  including  anything  more  than  adjacent  wavelengths 
in  the  space  charge  computation  is  that  they  can  be  expressed 
as  a geometric  profession  whose  'sum  to  infinity'  is  a sim- 
pler expression  (:j ■■_--)  than  the  sum  of  even  three  terms. 

The  assijmption  is  most  likely  to  become  unrealistic  at  the 
final  cavity,  where  the  next  wavelength  cihead  is  likely  to 
be  very  different  from  the  one  being  tracked,  if  the  effi- 
ciency is  high.  Ultimately  we  may  be  able  to  track  three 
consecutive  wavelengths,  the  outer  ones  acting  as  guards  for 
the  center  one. 


2.3  Subdivision  of  the  Beam 

Initially  the  one-wavelength  segment  of  the  beam  will  be 
divided  into  12  or  24  discs  and  these  will  be  tracked  for 
the  full  length  of  the  tube  to  establish  initial  values  of 
■Uie  rf  voltages  and  phases  at  each  gap,  using  an  existing 
disc  model  computer  program  [l].  We  shall  then  backtrack  to 
the  start  of  the  saturation  region,  subdivide  each  disc  into 
2,  3 or  4 concentric  rings,  and  repeat  the  calculation  from 
that  position;  v/ith  each  ring  now  moving  independently  under 
the  action  of  the  applied  fields.  Although  we  refer  to  these 
elements  of  the  beam  as  'rings',  we  do  not  think  of  them  as 
hydrodynamic  volume  elements  in  the  sense  that  Kosmahl  and 
Albers  [2]  consider  them.  In  this  work,  what  is  actually 
tracked  is  a 'super  electron'  having  about  10'  times  the 
charge  and  mass  of  a real  electron,  which  represents  the 
electrons  in  its  neighborhood.  Thus  discussion  of  'changes 
of  shape'  of  a ring  are  not  meaningful  in  this  context:  the 
ring  is  represented  by  a point  charge  which  has  no  shape,  but 
it  will  still  be  referred  to  as  a ring  for  brevity.  The  pre- 
cise charge  is  chosen  so  that,  when  multiplied  by  the  number 
of  rings  per  wavelength,  we  obtain  the  same  total  charge  as 
the  real  beam,  subject  to  a small  correction  to  be  discussed 
later. 


£23  'Three-Dimensional  Evaluation  of  Energy  Extraction  in 
Output  Cavities  of  Klystron  Amplifiers',  H.  G.  Kosmah]!. 
and  L.  U.  Albers,  IEEE  Transactions  on  Electron  Devices, 
ED-20  #10^  Oct.  1975,  pp.  883-890. 


2.4  The  Fields 


The  fields  acting  on  a ring  are: 

i)  the  rf  field; 

ii)  the  space  charge  field; 

iii)  the  magnetic  field. 

In  the  preliminary  disc  model  calcijlation  the  magnetic  field 
does  not  enter,  and  only  tl.e  axial  components  of  the  rf  and 
space  charge  fields  are  effective.  In  the  ring  model  part 
of  the  cal  dilation,  .ioth  axial  and  radial  components  of  all 
three  fields  are  to  Ije  included,  and  are  not  to  be  limited 
to  paraxial  approximations. 

It  will  be  noted  that  dc  electric  and  rf  magnetic  fields  are 
not  included;  the  effects  of  the  dc  electric  fields  in  the  gun 
are  represented  by  the  axial  injection  velocity  with  which 
the  electrons  are  started,  and  ’velocity- jump'  sections  are 
not  at  present  included.  Several  past  studies  have  shown 
that  the  rf  magnetic  fields  are  negligible  for  foreseeable 
microwave  tubes. 

There  are  various  methods  known  for  representing  the  fields 
in  computation.  They  may  be  derived  from  analytic  solutions 
of  the  wave  equation  or  Laplace’s  or  Poisson’s  equations  as 
appropriate,  or  from  Green’s  functions, or  from  the  gradients 
of  a potential  function.  We  have  available  a fast  trajectory 
algorithm  of  proven  accuracy  i_3j,  which  derives  the  fields  by 
interpolating  the  gradients  of  an  array  of  potentials  on  a 
rectangular  grid  which  overlays  the  interaction  region.  Thus 
our  working  representation  of  each  of  the  fields  will  be  a 
matrix  of  potentials  at  the  nodes  of  a suitable  grid.  TJie 
mesh  sizes  and  locations  of  the  grids  will  be  discussed  in 
detail  in  Section  2,6.  There  will  be  a separate  grid  and 
separate  matrix  for  each  of  the  three  fields. 


rjl  ’Electron  Ray- Tracing  Program  for  Image  Intensif iers ' , 
Final  Report,  Contract  DAAK02-67-C-0182,  by  J.R.M. 
Vaughan  and  0.  Buneman,  Sept.  1970. 
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2.5  Matrix  Representation  of  the  Fields 


For  each  field,  we  have  the  choice  of  constructing  either  a 
scalar  potential  matrix  or  a vector  potential  matrix;  for 
reasons  that  will  become  apparent  later,  we  choose  a matrix 
of  radius  x vector  potential  for  the  rf  fields,  a scalar 
potential  matrix  for  the  space  charge  fields,  and  a matrix 
of  vector  potential  x radius  for  the  magnetic  field.  These 
differences  are  not  apparent  to  the  ordinary  user,  but  must 
be  recognized  by  anyone  intending  to  delve  into  the  pro^am 
to  modify  it.  The  required  fields  (poxential  gradients}  are 
derived  from  a scalar  potential  matrix  by  differencing  the 
matrix  elements  in  the  same  direction,  but  from  a vector 
potential  matrix  by  differencing  in  the  perpendicular  direc- 
tion.  Thus  a scalar  (R,Z)  matrix  like  this: 


1 2 3 4 5 . . . 

12  3 4 5... 

1 2 3 4 5 . . . 

would  represent  a uniform  axial  field;  a ’vector  potential 
X radius ’ matrix  for  the  same  axial  field,  would  look  like 
this: 


0 0 0 

1 ' 1 1 

4 4 4 

9 9 9 


0 0 . . . 

1 1 . . . 

4 4 . . . 

9 9 . . > 


(In  practice,  of  course,  the  elements  are  not  simple  integers, 
and  the  scaling  factors  are  different  for  the  two  cases,  but 
the  vector  potential  matrices  dc  always  have  zeroes  along  the 
axis.)  It  may  be  worth  noting  here  another  possible  source 
of  confusion;  one  of  the  unfortunate  conventions  of  mathemat- 
ics is  that  matrices  are  printed  with  the  row  numbers  increas- 
ing downwards,  which  conflicts  with  Cartesian  coordinates  with 
y increasing  upwards.  Thus  we  shall  draw  meshes  superimposed 
on  the  interaction  space  of  the  tube  in  conventional  Cartesian 
form,  with  the  horizontal  lines  (representing  r rather  than  y) 
increasing  upwards.  But  in  a straight  printout  of  the  corre- 
sponding matrix,  the  top  line  of  the  matrix  will  correspond  to 
the  bottom  line  of  the  mesh,  and  vice  versa.  In  some  demon- 
stration cases  we  shall  deliberately  pi'ogram  the  computer  to 
print  a matrix  in  reverse  row  order  for  clarity,  but  a simple 
MAT  PRINT  statement  does  not  do  this. 
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If  it  is  later  decided  to  include  dc  electric  fields  to  rep- 
resent velocity  oiomp  sections,  a scalar  potential  matrix  will 
be  used  for  the  electrostatic  fields. 


2.6  Matrix  Dimensions 


The  fast  interpolation  routine  INTRA  for  the  potential  grad- 
ients requires  the  potentials  at  9 surrounding  mesh  points: 
thus  for  an  electron  at  Q in  Figure  3 , the  nearest  mesh  point 
is  P5,  and  the  remaining  points  P-j  to  P4  and  P5  to  P9  are  then 
determined  as  shown. 


Figure  3:  Mini-Matrix  for  INTRA 


The  routine  fits  an  exact  quadric  surface  through  these  9 
values,  and  obtains  the  gradients  of  the  two  principal  tan- 
gents at  Q to  th=  quadric,  representing  the  field  components 
at  Q.  (The  routine  INTRA  is  extremely  compact,  and  does  not 
explicitly  derive  the  quadric,  but  cuts  straight  through  to  the 
gradients,  without  neglecting  any  terms,  so  that  it  is  correct 
to  machine  accuracy.  It  was  derived  in  reference  [3] > where 
its  ad’.antage  over  5 point  interpolation  was  demonstrated.) 

One  can  see  from  Figure  3»  that  the  matrix  must  extend  at 
least  one-half  mesh  in  each  direction  beyond  any  position 
that  an  electron  Q can  occupy  during  the  calculation,  so 
that  9 surrounding  potentials  will  always  be  available. 

In  the  radial  direction,  an  electron  is  limited  by  the  tunnel 
wall  and  the  axis  (it  can  pass  through  the  axis,  but  its  rad- 
ial coordinate  is  by  definition  always  positive,  so  that  it 
appears  in  the  R-Z  plane  to  bovuice  off  the  axis) . Thus  the 


minimum  radial  mesh  system  would  extend  from  one-half  mesh 
below  the  axis  to  one-half  mesh  above  the  wall.  But  there 
is  such  obvious  convenience  in  having  one  of  the  mesh  lines 
along  the  axis,  and  another  along  the  wall,  that  we  choose 
to  make  every  matrix  (for  rf,  space  charge  and  magnetic 
fields)  extend  fadially  from  1 inesh  below  the  axis  to  1 mesh 
above  the  wall,  recognizing  that  this  makes  the  radial  matrix 
dimension  greater  by  1 than  it  would  strictly  have  to  be. 

The  number  of  meshes  between  the  axis  and  the  tunnel  wall 
need  not  be  the  same  for  all  three  matrices.  For  the  and 
magnetic  field  matrices,  the  mimbers  may  be  chosen  at  vill  — 
the  larger  the  number,  the  more  accurate  can  be  the  represen- 
tation of  the  field,  but  the  larger  is  the  memory  requirement, 
and  the  more  computation  is  required  to  set  up  the  matrix. 

It  does  not,  however,  affect  tha  amount  of  computation  in  the 
main  ring- tracking  part  of  the  program:  at  each  step  9 adja- 
cent values  have  to  be  extracted  and  interpolated,  and  it 
makes  no  difference  whether  they  are  9 out  of  100  or  9 out 
of  1000.  Typical  values  for  the  number  of  radial  meshes  will 
range  from  4 for  rough  calculations  or  debugging,  to  about  20 
for  precise  work  (there  is  no  real  advantage  in  going  to  rad- 
ial mesh  numbers  that  are  higher  than  the  ratio  of  tunnel 
radius  to  ferrule  comer  radius,  which  is  typically  not  more 
than  about  20).  The  radial  mesh  numbers  are  denoted  NgR  for 
the  rf  vector  potential  and  Nwm  for  the  magnetic  vector  poten- 
tial matrices.  Allowing  for  ‘tne  guard  rows,  the  matrices  run 
from  -1  to  Ngo  + 1,  and  -1  to  + 1.  When  the  program  calls 
for  ’mesh  numbers',  it  is  the  basic  numbers  Nqr,  Nf/iR,  etc. 
that  are  to  be  entered,  ''he  program  will  add  the  guard  rows 
and  columns  as  necessary. 

For  the  space  charge  matrix,  the  fast  algorithm  to  be  given 
in  Section  5 requires  that  the  number  NgR  of  radial  meshes 
be  a power  of'  2.  It  will  usually  be  4 or  8,  possibly  16.  In 
this  case  the  choice  does  affect  the  main  computation  speed, 
since  this  entire  matrix  has  to  be  recalculated  after  every 
time  step,  and  this  is  the  pacing  item  for  the  whole  program. 

In  the  axial  direction,  the  number  of  meshes  is  similarly  a 
free  choice  for  the  rf  and  magnetic  field  matrices  (except 
that  it  must  be  an  even  number  for  the  rf  matrix;  an  odd  num- 
ber would  be  an  unlikely  choice  for  either  matrix  in  any  case) 
But  the  choice  is  strictly  limited  to  6,  12,  24  or  possibly  48 
for  the  space  charge  matrix.  The  rf  mesh  is  physically  tied 
to  the  cavity  period,  as  shown  in  Figure  4(b),  but  with  an 
extra  mesh  at  each  end.  This  matrix  is  stationary,  but  is 


repeated  for  every  cavity.  The  matrix  represents  the  poten- 
tials for  1 volt  peak  rf  across  the  gap  at  zero  phase,  and 
in  use  the  gradients  will  be  multiplied  by  appropriate  volt- 
age and  phase  factors  for  each  cavity.  If  the  cavity  period 
is  divided  into  Nqa  parts,  the  matrix  numbering  will  run  from 
-1  to  Nca  ^ 1*  Nca  will  typically  be  not  less  than  4 nor 
more  than  50. 

The  magnetic  field  matrix.  Figure  4(c),  will  similarly  be 
tied  to  the  magnet  period.  For  ’single  period'  focusing 
(alternating  magnet  polarities  in  each  cavity)  this  is  twice 
the  cavity  period;  for  'double  period'  focusing,  the  magnet 
period  is  four  times  the  cavity  period.  There  is  also  the 
possibility  of  focusing  systems  being  used  which  do  not  tie 
the  magnet  period  to  the  cavity  period  at  all,  so  we  shall 
allow  the  magnet  period  to  be  an  independent  variable,  but 
with  the  expectation  that  in  most  cases  it  will  be  specified 
as  2 or  4 times  the  cavity  period.  An  example  of  a nonuni- 
form field  that  was  not  tied  to  the  cavity  period  would  be 
a field  produced  by  a solenoid  of  several  coils  with  inde- 
pendent current  controls,  so  that  a ' programmed ' field  could 
be  generated.  This  would  be  treated  as  a periodic  field 
whose  period  extended  over  all  the  saturation  region  cavities, 
so  that  the  computation  would  never  get  beyond  the  first  per- 
iod. If  the  magnet  period  is  divided  into  Nm  parts,  the 
magnet  matrix  will  run  from  - 1 to  Nma  + 1 , and  will  be  re- 
peated for  every  magnet  period.  Typical  values  of  will 
be  from  4 to  24.  For  both  rf  and  magnetic  field  matrices, 
there  is  a two  mesh  overlap  of  consecutive  matrices,  but 
there  is  no  confusion  as  to  which  one  is  to  be  used  for 
rings  in  the  overlap  range;  if  a ring  is  on  or  to  the  right 
of  the  tunnel  midplane,  it  uses  the  matrix  on  the  right;  if 
it  is  to  the  left  of  the  midplane,  it  uses  the  matrix  on  the 
left. 

If  the  magnetic  field  is  uniform  (solenoid  focusing),  it  is 
not  necessary  to  construct  a magnetic  matrix  at  all;  the 
trajectory  program  will  allow  for  a uniform  field  by  analytic 
methods . 

The  space  charge  matrix  is  different  in  character:  the  number 
of  meshes  in  one  beam  wavelength  must  be  one  of  the  nxambers 
for  which  a superfast  F.F.T.  exists;  usually  it  will  be  12 
or  24;  and  the  corresponding  grid  is  not  stationary  but  is 
moving  with  the  beam.  (Fig.  4(d))  If  the  number  of  meshes 
is  NsAf  we  can  arrange  that  the  mesh  position  of  a ring  is 
always  within  the  range  0.5  to  N5;^  + .5  (since  we  have  al- 
ready agreed  that  a ring  can  be  moved  up  or  back  1 beam  wave- 


length  to  keep  it  in  range),  so  that  this  matrix  runs  from 
0 to  Ng^  + 1 . 

In  the  Super  Basic  debugging  version  of  the  pro^am,  the 
matrices  can  be  dimensioned  exactly  as  written  (e.g.,  -1 
to  9 radially.  0 to  13  axially  for  a nominal  8x12  space 
charge  matrix),  but  for  the  working  Fortran  version  there 
is  the  added  complication  that  zero  or  negative  indices  are 
inadmissible,  so  that  all  the  indices  have  to  be  shifted  up- 
ward by  1 or  2 as  the  case  may  be.  This  is  a thorough  nui- 
sance, and  it  is  to  be  hoped  that  eventually  a version  of 
Fortran  will  come  out  that,  like  PL/1  and  Super  Basic,  allows 
negative  or  zero  indices.  In  the  meanwhile,  it  is  a quirk  to 
be  recognized  by  anyone  digging  into  the  program  details,  but 
irrelevant  to  the  ordinary  user. 

Since  the  axial  and  radial  mesh  sizes  for  each  matrix  are 
determined  independently,  the  meshes  will  not  in  general  be 
square.  Normally  the  choices  made  are  such  that  they  are 
elongated  in  the  axial  direction;  the  interpolation  routine 
allows  for  this,  but  there  is  some  advantage  in  making 
choices  that  do  not  result  in  extreme  elongation  — say, 
not  more  than  8:1.  It  does  not  appear  likely  that  a case 
would  ever  arise  in  which  the  meshes  were  elongated  radially. 


2.7  Coordinate  Systems 

The  basic  coordinate  system  of  the  program  is  a stationary 
Cartesian  system  in  MKS  units:  the  Z axis  lies  along  the 
tube  axis,  and  the  origin  is  at  the  tunnel  midplane  on  the 
entrance  side  of  cavity  #1 . The  Z coordinate  of  a disc  or 
ring  at  any  time  is  its  distance  in  meters  from  this  plane, 
and  the  R coordinate  of  a ring  is  similarly  in  meters.  R 
will  be  broken  down  into  X and  Y components  in  the  trajectory 
computation,  the  XZ  plane  being  the  plane  initially  contain- 
ing a super-electron.  (It  will  move  out  of  this  plane  unless 
the  magnetic  field  is  zero  everywhere.)  The  XYZ  axes  remain 
fixed,  but  each  super-electron  has  its  own  RZ  plane  which 
rotates  about  the  Z axis  so  that  it  always  passes  through 
the  current  position  of  the  super-electron,  as  shown  in 
Figure  5. 

Each  of  the  three  matrices  constitutes  an  aioxiliary  coordi- 
nate system  in  which  the  units  are  the  mesh  sizes.  For  each 
ring,  we  shall  know  from  its  Z coordinate  which  cavity  and 
which  magnet  period  it  is  in,  so  we  shall  subtract  the  Z 
coordinate  of  the  origin  of  the  matrix  for  that  cavity  or 
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Figure  5:  Mesh  line  numbering  and  Node  numbering 

for  INTRA.  Three  different  meshes  are 
used  for  the  rf  vector  potential,  the 
magnetic  vector  potential,  and  the  space 
charge  potential,  but  all  are  numbered  in 
the  same  way. 
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magnet,  and  divide  by  the  axial  mesh  size  to  get  the  relative 
position  in  mesh  units. 

Similarly  'the  radial  position  in  mesh  units  is  the  MKS  r 
divided  by  the  radial  mesh  size,  usually  different  from  the 
axial  mesh.  The  relative  position  in  mesh  imits  then  allows 
us  to  obtain  the  gradients  representing  that  iield. 

The  same  procedure  applies  to  the  moving  grid  in  which  the 
space  charge  forces  are  evaluated.  This  grid  moves  with  the 
'dc  beam  velociti-'  , but  this  terra  is  slightly  ambiguous  when 
potential  depression  is  taken  into  account.  For  reasons 
given  in  Hj  we  choose  the  velocity  at  a radius  b/'/3  in  the 
initial  uniform  beam  as  the  nominal  dc  beam  velocity.  The 
moving  grid  is  assigned  this  velocity,  and  retains  it  through- 
out the  motion,  so  that  in  saturation  the  beam  is  mostly  slid- 
ing back  through  it.  The  zero  of  time  is  the  instant  at  which 
the  origin  of  the  moving  grid  passes  the  origin  of  the  fixed 
MKS  coordinate  system.  Since  the  zero  of  the  moving  grid  is 
at  its  left-hand  end,  this  implies  that  the  beam  segment  to 
be  tracked  crossed  the, entrance  plane  (mid- tunnel  on  the  left 
of  cavity  #1)  before  t=0,  and  is  already  distributed  through 
cavity  1 and  part  of  cavity  2 att=0  (the  1 wavelength  beam 
segment  is  typically  about  1-3A  cavities  long). 

It  is  evident  that  this  use  of  four  separate  coordinate  sys- 
tems involves  an  enormous  number  of  transformations,  but  they 
are  extremely  simple  and  fast  operations  on  the  computer.  To 
compel  all  the  fields  to  use  a common  set  of  mesh  units  would 
force  undesirable  compromises  on  all  of  them.  By  letting  each 
grid  be  independent,  and  detemiined  only  by  its  own  constraints 
and  accuracy  needs,  while  relating  each  to  the  underlying  MKS 
coordinates,  we  retain  great  flexibility,  and  freedom  to  in- 
corporate additional  matrices,  such  as  one  for  electrostatic 
fields  if  we  want  to.  To  this  writer  at  least,  there  is  also 
a strong  psychological  advantage  in  using  I/KS  units  as  the 
basic  system,  rather  than  normalized  units  such  as  Pierce's 
y>  it  gives  a feeling  of  knowing  where  the  ele'^trons  ' really 
art ' . Certainly  if  programming  errors  or  incorrect  data  en- 
tries result  in  unrealistic  values,  this  becomes  much  more 
obvious  if  they  are  expressed  in  familiar  units. 
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2.8  Input  Dimensions 

Since  tube  drawings  are  still  often  dimensioned  in  inches, 
the  user  is  allowed  to  make  a choice  of  entering  all  linear 
dimensions  in  either  inches  or  millimeters.  The  choice  sets 
a conversion  factor  CLIN  to  either  .0254  or  .001  respective- 
ly. Once  made,  the  choice  must  be  adhered  to  for  all  inputs 
involving  linear  measures . 

The  program  then  converts  all  lengths  and  distances  to  meters 
by  multiplying  by  CLIN.  Output  coordinates  are  converted  to 
millimeters,  but  conversion  back  to  the  input  units  could  be 
substituted  very  easily  if  this  is  preferred.  This  use  of 
the  most  familiar  \mits  for  input  is  considered  of  great 
importance  for  avoiding  wasted  runs  caused  by  incorrect  data 
entries. 

An  example  of  the  very  straightforward  input  for  the  prelim- 
inary time- sharing  version  of  the  program  is  shown  in  Figure 
6.  The  user  needs  to  know  the  physical  parameters  of  the  tube 
he  wishes  to  have  calculated,  and  to  have  some  idea  of  the 
accuracy  level  he  wants,  to  allow  a suitable  choice  of  matrix 
dimensions,  but  he  need  know  nothing  more  about  the  program. 

The  FORTFIAN  input  is,  as  always,  somewhat  more  restricted  in 
format,  but  is  fully  described  in  the  User's  Manual. 
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Figure  6:  Program  Input 
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3.0  MAGNETIC  FIELD 


The  only  restriction  placed  on  the  magnetic  field  is  that  it 
be  axisymmetric,  that  is,  that  it  have  no  dependence  on  0. 

In  PPM  structures  the  field  becomes  purely  radial  at  certain 
planes,  and  of  course  j.t  is  purely  axial  at  the  centers  of 
the  gaps;  if  we  are  to  model  this  complete  range  of  directions 
accurately,  we  cannot  allow  any  paraxial  approximations.  The 
method  to  be  described  passes  from  purely  axial  to  purely 
radial  with  no  loss  of  accuracy,  and  constitutes  a valid  solu- 
tion of  Laplace’s  equation.  (Some  published  approaches  to  this 
problem  use  ad  hoc  expressions  v/hich  do  not  satisfy  Laplace  f^j)* 

Primarily  the  magnetic  field  is  represented  in  the  computer 
program  by  the  parameters  of  a set  of  ideal  circular  current 
loops,  usually  not  more  than  10  in  number.  These  loops  are 
chosen  so  that  the  field  they  generate  matches  the  actual  field 
over  the  working  region  within  a desired  tolerance.  However, 
as  explained  in  Section  2.5,  the  trajectory  algorithm  requires 
d matrix  whose  elements  are  ’radius  X magnetic  vector  potential' 
at  the  nodes  of  a suitable  mesh  in  the  r-z  plane.  Therefore, 
the  loop  parameters  are  used  to  generate  this  matrix,  which  then 
becomes  the  working  representation  of  the  field  for  the  main 
ray- tracing  part  of  the  program.  Whether  the  chosen  loops  rep- 
resent the  desired  field  accurately  or  not,  the  field  derived 
from  the  loops  is  always  a true  solution  of  Laplace’s  equation. 

The  reasons  for  this  choice  of  method,  and  its  implementation, 
will  now  be  discussed  in  more  detail. 


3.1  The  Vector  Potential  Matrix 

The  ray-tracing  x'outine  derives  information  about  the  magnetic 
field  by  extracting  the  values  at  the  9 nearest  nodes  of  a 
potential  matrix  for  each  ring  at  each  time  step.  These  values 
are  then  interpolated  by  subroutine  INTRA  to  give  the  gradients 
at  the  position  of  the  ring.  Since  the  magnetic  field  is  static, 
a scalar  potential  matrix  could  be  used.  But  it  is  found  that 
the  magnetic  field  terms  can  be  more  effectively  integrated  into 
the  ray-tracing  routine  if  the  vector-potential  is  used,  with 
each  value  multiplied  by  R.  The  vector  potential  is  a less 


£43  H.K.  Detweiler  and  J.E.  Rowe,  'Electron  Dynamics  and 
Energy  Conversion  in  0-Type  Linear  Beam  Devices’  in 
'Advances  in  Microwaves’,  Vol.  / , 1971 » Academic  Press, 
p.  35.  The  pair  of  equations  (14)  on  p.  39  do  not 
satisfy  Laplace. 
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familiar  concept  but  it  takes  a very  simple  form  for  axi- 
symmetric  fields,  and  is  easily  calculated  from  formulas  to 
be  given.  It  should  be  remembered  that  the  use  of  vector 
potential  involves  vector  cross  products,  so  that  the  radial 
differences  in  the  matrix  determine  the  axial  field  and  vice 
versa. 


3.2  Tlie  Magnetic  Vector  Potential 


For  a general  axisymmetric  field,  the  vector  potential  at 
(R,Z)  is 


R 


r dr 


= ^ / 
o 

= • (flux  through  circle  of  Radius  R) 


The  direction  of  the  vector  A is  circumferential.  The  quan- 
tity to  be  stored  in  the  matrix  element  corresponding  to 
(r,z)  is  M = rA,  and  the  field  components  are  then  given  by 


B , B = 

z r dr  ’ r 


1 M 

r dz 


3.2.1  Uniform  Field 


If  the  field  is  uniform 


B = B 
z o 


B = 0 
r 


independent  of  r 
everywhere . 


Then 


M(R,Z)  = R.A  = R^Bq/2 


(3-1) 


(3-1a) 


(3-2) 


Thus  the  matrix  for  a uniform  field  region  is  easily  recog- 
nizable: in  each  column,  the  values  are  proportional  to  0,  1, 
4,  9 I etc.,  and  all  columns  are  identical. 
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However,  for  the  uniform  field  case,  the  ray- tracing  program 
will  be  diverted  to  a simpler  set  of  equations  in  which  only 
the  axial  component  B2  occurs,  so  that  it  is  not  necessary  to 
construct  the  M matrix  at  all.  This  includes  the  case  of  zero 
field. 


3 • 2 . 2 General  Axisymmetric  Field 

The  basic  sources  of  an  axisymmetric  magnetic  field  are  loops 
of  current  flowing  in  planes  perpendicular  to  the  axis.  If 
the  field  is  generated  by  a solenoid,  there  is  aoi  obvious  cor- 
respondence between  these  mathematical  current  loops  and  the 
actual  turns  of  the  coil.  But  if  the  physical  source  is  a 
permanent  magnet,  one  can  still  visualize  the  billions  of 
'Amperian'  currents  circulating  in  the  aligned  molecules,  can- 
celling each  other  everywhere  in  the  interior  (for  uniform 
magnetization),  but  adding  up  to  a large  surface  current  dens- 
ity at  free  surfaces  that  are  not  pei'pendicular  to  the  direc- 
tion of  magnetization  (Figure  7).  It  should  be  noted  that  the 
kind  of  ring  magnet  often  used  in  TWTs  has  two  such  surfaces, 
the  inner  and  outer  cylindrical  surfaces,  with  Amperian  cur- 
rents in  opposite  directions.  The  correct  representation  of 
this  magnet  therefore  requires  two  sets  of  current  loops  of 
opposite  polarity  located  on  the  inner  and  outer  diameters. 

The  writer  has  seen  quite  large-scale  attempts  to  compute 
fields  based  on  the  assumption  that  the  field  can  be  repre- 
sented by  a single  current  sheet  or  set  of  loops.  This  is 
only  true  if  one  confines  attention  to  the  region  close  to  the 
axis:  we  shall  find  that  in  this  case  (which  is  common  in  TWTs 
of  course)  a single  set  of  coils  can  be  sufficient;  but  it 
should  be  remembered  that  this  is  not  generally  true. 

The  most  common  textbook  expression  for  the  field  is  an  ex- 
pansion in  terms  of  the  axial  values  and  their  differentials. 
Evidently  the  writers  of  textbooks  have  never  actually  carried 
out  this  calculation  because,  while  algebraically  sound,  the 
method  has  numerical  instabilities  which  make  it  useless  in 
practical  cases  ^5].  It  only  works  if  one  restricts  oneself 
to  paraxial  cases, "a  limitation  which  we  have  specifically  re- 
jected, or  if  the  axial  values  of  the  field  are  known  with 
machine  accuracy  (order  of  8 digits  or  more).  The  values  can 
only  be  known  with  this  kind  of  accuracy  if  one  has  derived 
them  mathematically  from  the  'sources'  just  described,  which 
implies  that  one  knows  what  the  sources  are.  If  one  does  know 


[3l  J.R.M.  Vaughan,  'Representation  of  Axisymmetric  Magnetic 
^ Fields  in  Computer  Programs',  IEEE  Transactions  on  Elec- 

tron Devices,  ED-19  #2,  February  1972,  pp.  144-151. 
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Figure  7 

Permanent  Magnet  and  Equivalent  Surface  Currents 


this,  then  it  is  much  more  natural  to  derive  the  desired 
matrix  M directly  from  the  sources  than  via  the  axial  dif- 
ferentials. If  what  one  knows  is  really  a set  of  data  values 
of  the  field,  then  the  procediire  will  be,  first  to  find  a set 
of  sources  — ideal  current  loops  --  which  will  represent  the 
data,  and  then  to  derive  M as  before.  Several  methods  of  find- 
ing appropriate  loop  parameters  in  practical  cases  have  been 
described  [6]. 

Since  the  field  is  axi symmetric,  every  current  loop  is  cen- 
tered on  the  axis  and  perpendicular  to  it;  each  has  three 
parameters:  axial  position  Z^;,  radius  and  current  the 
writer  prefers  to  use  the  'strength'  M-]  = jigl'i/Zn  as  an 
equivalent  parameter. 

The  axial  positions  of  the  loops  are  unrestricted:  the  grid 
covers  one  period  of  the  magnet  structure  in  the  axial  direc- 
tion, but  loops  lying  outside  that  axial  range  can  be  con- 
tributing to  the  field  within  the  range.  For  'single  period' 
focusing,  for  example,  the  period  covers  two  consecutive  cavi- 
ties; the  field  is  represented  by  four  coils,  two  in  the  gaps 
of  these  cavities,  and  one  in  the  next  cavity  gap  on  either 
side,  so  the  two  latter  have  z positions  outside  the  range  of 
tlie  grid. 

The  radius  of  a coil  restricted:  physically,  it  must  be 
greater  than  the  tunnel  radius,  otherwise  an  electron  could 
encounter  a field  singularity.  This  condition  is  in  practice 
only  violated  if  a mistake  has  been  made  in  calculating  (or 
entering)  the  coil  data.  But  when  we  use  a potential  grid  to 
represent  the  field,  a somewhat  stronger  condition  is  required: 
the  coil  must  lie  not  only  outside  the  turmel,  but  beyond  the 
outermost  grid  line  by  about  0.5  mesh  so  that  no  mesh  point 
can  lie  too  close  to  the  singularity.  In  practice,  the  correct 
coil  position  for  a typical  PPM  struocure  is  at  about  1.4  or 
1.5  times  the  tunnel  radius,  so  this  condition  only  comes  into 
play  if  a very  coarse  mesh  is  used  = 2 or  3,  for  example j. 

The  program  checks  each  coil  radius,  and  if  it  is  too  small  for 
the  mesh  size  chosen,  it  will  automatically  increase  the  radius 
by  a factor  C5  to  bring  it  to  the  minimum  acceptable  value.  It 
simultaneously  increases  the  strength  by  a factor  C5(1+.25C5(C5-1 ) ) , 
which  restores  the  field  strength  on  the  axis  to  the  correct 
value.  A diagnostic  is  printed  specifving  the  new  values  assigned. 


[6j  J.R.M.  Vaughan,  'Methods  of  Finding  the  Parameters  of 
Ideal  Current  Loops  for  Computer  Simulation  of  Magnetic 
Fields',  IFee  Transactions  on  Electron  Devices,  ED-21 
#5,  May  1974,  pp.  310-312. 
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For  a nonuniform  solenoid  focused  case,  the  coil  radii  will 
all  be  much  larger,  and  the  problem  would  not  arise  even  at 

N„r  = 2. 

In  the  'inch-gauss'  or  'millimeter- gauss'  units  to  be  used 
for  data  input  (‘see  Section  3.3),  the  'strength'  of  a single 
loop  of  radius  inches  or  mm  generating  a field  of  G-j  gauss 
at  its  center  is  R^G^/n. 


The  flux  through  the  circle  R,Z  due  to  the  source  R^ , 
is  [7j  111 


$ = UqI^(RR.|)  = 


£ - C K(c)  - ^ E(c) 


where  c^  = 4RR^  /UR+R^)^  + .(Z-Z^)' 


and  K and  E are  the  complete  elliptic  integrals,  modulus  c. 
(The  alternative  expansion  in  Legendre  polynomials  has  very 
slov;  convergence  over  much  of  the  range  we  shall  need.) 


Combining  (3-3)  with  (3- 1a)  to  obtain  A,  and  multiplying  by 
R we  have 


M(R,Z) 


^0^1 

~7T 


(RRl  )^ 


c)  K(c)  - I E(c) 


where  the  reason  for  using  p I^/2tt  as  a coil  parameter  is 
now  evident. 


For  values  of  c s 0.2,  (3-5)  can  be  evaluated  by  using  the 
elliptic  integral  subroutine  ELIVA,  which  is  incoi’porated 
in  the  program.  It  is  significantly  faster  than  the  IBM 
routine.  For  c <0.2,  the  terms  in  the  square  bracket  be- 
come nearly  equal,  and  we  improve  the  accuracy  by  replacing 
them  with  the  power  series  expansion 


TTC^ 

IF 


3 2 


4 245  6] 

c + ^ c 


r?]  J.  Jeans,  'The  Mathematical  Theory  of  Electricity  and 
Magnetism' , Cambridge  Univ.  Press,  5th  Ed.  1933,  p.  443. 
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The  last  term  is  <10  for  c < .2,  so  further,  teims  are  \m- 
necessary. 

3.3  Prof;ram  Input 

In  line  with  the  policy  of  inputting  data  in  familiar  or  con- 
venient units,  the  coil  data  will  be  called  for  in  inch-gauss 
or  mil.. imeter- gauss  units.  These  are  then  converted  to  MKS 
units  vs.lng  CLIN  = .0254  or  .001,  ai\d  CMAG  = .0001  (for  gauss 
to  Tesla).  Note  that  (3-4)  cannot  be  interpreted  directly  in 
mesh  units,  because  the  mesh  units  normally  differ  for  R and 
Z;  it  could,  of  course,  be  modified  to  allow  for  this,  but  it 
appears  simpler  to  keep  the  R's  and  Z‘s  in  meters  and  evaluate 
(3-4)  and  (3-5)  as  written. 

Thus,  if  the  magnet  period  is  inches,  and  is  di-vided  into 
N^.  meshes  axially,  the.Z  values  are 


* CLIN  * I/N^ 


for  I = 1 to  N, 


or  for  I = 1 to  Nj^j^/2 


if  the  magnet  period  has  Z symmetry. 

Similarly  if  the  tunnel  radius  is  Aj  inches,  and  is  divided 
into  Nj^j^  meshes,  the  R values  are 


Aj  * CLIN  * J/Nj 


for  J = 1 to  Nj^  + 1 . 


(3-8) 


Subroutine  MAGVA  evaluates  (3-5),  using  (3-4)  and  (3-6) 
where  appropriate,  for  these  ranges  of  I and  J,  for  each 
coil.  The  M values  for  several  coils  are  additive:  although 
they  are  strictly  vectors,  they  all  have  the  same  direction 
(normal  to  the  R-Z  plane)  so  they  can  be  added  algebraically. 
The  PqIi/2tt  terms  in  (3-5)  are  the  entered  strengths  M,  mul- 
tiplied by  CMAG  x CLIN. 

The  number  of  coils  necessary  to  represent  one  period  of  the 
field  has  never  so  far  exceeded  10;  4 coils  are  sufficient 
for  ordinary  'single  period'  focusing,  and  8 for  double  period. 
In  either  case,  only  half  the  matrix  need  be  calculated,  the 
other  half  being  the  same  with  reversed  signs.  For  cases 
such  as  the  multi-coil  nonuniform  solenoid,  the  whole  matrix 
must  be  calculated. 
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The  remaining  elements  of  the  full  matrix  can  then  be  filled 
in  without  further  calculation: 


M (R, 

0)  = 

M (R, 

^MA^ 

for 

R = 

1 

to 

M(R, 

-1)  = 

M (R, 

Nma-1) 

for 

R = 

1 

to 

”mr  * 

M(R, 

) = M 

(R,  1) 

for 

R = 

1 

ta 

®*MR  * '' 

M (0, 

z)  = 

: 0 

for 

Z = 

- 

1 to 

«MA  * ' 

M (-1 

, Z) 

= M (1 

, z) 

for 

Z = 

- 

1 to 

»«A  * 1 

The  program  will  allow  three  options  for  the  magnetic  poten- 
tial matrix: 


i)  Compute  the  matrix  and  discard.it  at  the  end  of  the 
run. 

ii)  Compute  the  matrix  and  save  it  on  a file  MAGMAT  for 
future  use . 

iii)  Read  in  the  matrix  from  MAGMAT. 


The  nominal  matrix  dimensions  and  Njyj.  are  stored  with  the 
matrix,  and  in  case  (iii)  are  compared  with  the  specified 
values  as  a safeguard  against  reading  in  an  incorrect  matrix. 

The  complete  process  of  generating  the  magnetic  vector  poten- 
tial matrix  from  the  original  specification  of  the  field  is 
summarized  in  the  flow  chart  in  Figure  6. 

There  remains  the  question  of  location  of  the  magnet  period 
in  relation  to  the  cavities:  the  matrix  is  only  needed  for  the 

last  k cavities  in  which  the  ring  model  of  the  beam  is  to  be 

used  (k  10).  The  convention  adopted  is  that  the  Nwj^'th 
grid  line  of  the  last  magnet  period  coincides  with  the  mid- 
plane in  the  tunnel  following  the  last  cavity.  Then  as  many 
repetitions  of  the  magnet  period  are  added  prior  to  this  as 
are' necessary  to  extend  back  over  at  least  k+1  cavities  (since 
if  we  change  from  the  disc  to  the  ring  model  at  cavity  k,  some 
elements  of  the  beam  segment  will  still  be  back  in  cavity  k-1). 

For  example,  if  the  tube  has  50  cavities,  of  which  the  last  10 

are  to  be  computed  with  the  ring  model,  and  if  'double  period' 
focusing  is  used  (magnet  period  = 4 cavity  periods),  then  mag- 
net period  1 will  cover  cavities  39-42,  period  2 cavities  43-46, 
and  period  3 cavities  47-50.  If  12  cavities  were  to  be  used, 
then  4 magnet  periods  would  be  needed,  starting  at  cavity  35. 
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A:  alternative  starting  points,  depending  on  how  the  magnetii 
field  is  specified. 


Figure  8.  Flow  chart  for  Magnetic  Vector  Potential  computation. 
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3.4  Numerical  Example 

Since  it  is  rather  easy  to  become  confused  over  the  units, 
and  is  useful  to  have  a test  case  for  debugging,  a fully 
worked-out  numerical  example  follows,  based  on  the  standard 
test  case  of  Sec'tion  2.1.  The  pertinent  data  are: 


Tunnel  diameter 

.200" 

FerrifLe  O.D. 

. 300" 

Magnet  I .D. , O.D. 

1.5", 

Double  cavity  length 
(i.e.,  magnet  period) 

. 600" 

Magnet  length 

.200" 

Ferrifle  gap 

.1" 

Web  thickness 

• .1" 

Magnet  material : Samarium 

Cobalt 

This  information  is  sufficient  for  calculation  of  the  field 
by  the  method  of  Sterzer  and  Siekanowicz  [8]  which  is  em- 
bodied in  Litton  proprietary  program  /PPMMAG18/.  Figure  9 
shows  the  case  run  on  this  program.  Half  way  down  the  page 
we  see  the  gap  center  field  midway  betv/een  the  ferrules, 

5316.9  gauss,  and  the  field  on  the  axis  in  the  same  Z plane, 
3046.7  gauss.  These  two  values  are  sufficient  for  a first 
approximation  to  the  parameters  of  a coil  to  produce  this 
field,  using  the  method  of  [5j  and  :‘61  . The  program  then 
assigns  identical  coils  of  opposite  polarity  in  adjacent  gaps, 
and  a fourth  coil  beyond  the  third,  and  calculates  the  pertur- 
bations they  cause  in  the  first  gap. ' The  coil  radii  and 
strengths  are  then  adjusted  to  restore  the  fields  at  r=0  and 
r=a  to  the  desired  values,  and  the  two  outer  coils  are  slight- 
ly adjusted  to  force  the  field  to  zero  at  the  mid- tunnel  posi- 
tions; the  final  parametei-s  for  all  four  coils  are  printed  out 
at  the  bottom  of  the  page  in  the  desired  inch-gauss  units. 

The  program  gives  other  information  which  is  useful  for  magnet 
design  but  not  relevant  to  the  present  application. 


[8]  F.  Sterzer  and  W.W.  Siekanowicz,  'The  Design  of  Periodic 
Permanent  Magnets  for  Focusing  of  Electron  Beams',  RCA 
Review,  Vol.  18,  pp.  39-59,  Mar.  1957. 
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MAGNET  OPERATING  POINT.  H=-F.G58  K-OE.  E=  4.913  KG»  E/H=  -1.349 


1. 

2. 

3. 

4. 

G. 

5. 


REGION  1 

PERMEANCE  <UL> 

FLUX  <Gl 

TUNNEL 

.261 

894 

FERRULE-FERRULE 

.393 

1347 

FERRULE  O.D. 

.309 

1 061 

FERRULE  TO  STEP 

8.448 

28971 

STEP  TO  MAGNET 

. 000 

0 

EXTERNAL 

3. .290 

1 1282 

TOTAL 

12.700 

43555 

GAP  CENTER  FIELD  = 531 G. 9 GAUSS 
AXIAL  PEAKS  AND  MINIMA: 

Z GAUSS 

.1500  304G.7 

.3000  .0 

R.M.S.  FIELD  1873.2  GAUSS 


INTEGRATED  FIELD  <Y.N)  ? H 


HAPNDNIC  AMPLITUDES: 

MO.  AMPLITUDE  PCT.  OF  PEAK 
1 £816.8  85.9 

3 -411.9  -13.5 

5 20.1  .7 

7 2.1  .1 

FLUX  DEMSITY  IM  SHIM  AT  D3  21236  GAUSS 

AT  D2  1 0363 
IM  FERRULE  13035 


WEIGHT  PER  RF  CAVITY  = .1627  LBS 

WEIGHT  OF  ONE  MAGNET  = .0966  LBS  < 43.837  GRAMS' 


CALCULATE  E'X 'VALENT  COILS  <Y.H)  ? Y 


GAP/'AXIS  FIELD  PATIO  1.745*  COIL  FAD  PATIO  1.393 
FIRST  VALUE;-  COIL  RAD  .139.  STRENGTH  135.1 
SECDHD  ESTIMATE  OF  RADIUS  .14801348 


PMS  DF  COIL  FIELDS  2010.5  GRiJSS 
AXIAL  DATA  MOW  ON  FILE  "VGIPLOT/ 


COIL  DATA  IN  INCH-GAUSS  <1>  OP  MESH  <M'  UNITS  ? I 


COIL  NO.  1 

rrm  nn,  ^ 

COIL  no!  3 
COIL  NO.  4 


RADIUS 
. 1480 

, 1 J80 
! 1480 
.1480 


AXIAL  rns. 

-.1500 

. J50rt 

I45OO 

.7500 


STRENGTH 

-158.2 

169.8 

-169.8 

158.2 


AtnVE  VALUES  TO  PE  DIVIDED  tV  'UNITIN'  FOP  U:E  IN  SLAC  PROGRAM 


I 


INPUT  1 FOP  NEW  value:.  ? FOP  A NFM  PUN.  OR  3 TO  STOP 
IF  NFW  VAIUFS.  INPUT  THFM.  THEN  TYPE  'GO'  3 r 
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This  set  of  coils  reproduces  the  desired  field  values  within 
0.^%  at  r=0  and  .3%  at  r=a,  at  mid-gap.  The  rms  value  of  the 
coil  fields  over  a complete  period  is  about  7%  high.  If  nec- 
essary this  discrepancy  could  be  reduced  by  using  more  coils 
to  represent  th'-  field,  but  this  does  ' seem  necessary. 
Other  methods  of  obtaining  coil  para*-  ters  without  using 
/PPMMAG18/  are  detailed  in  {_3J  an''  6J  . 

The  input  data  relevant  to  the  magnetic  vector  potential  ma- 
trix is,  for  our  standard  test  case: 


Tunnel  diameter  .2" 

Magnet  period  .6”  ^M  ~ 

Radial  mo  s hes  in  tunnel  radius  = 4 

Axial  meshes  in  magnet  period  Nj^  = 16 
Coil  data  R1 , Z1  and  M1  exactly  as  in  Figure  9. 


Lm  = .6" 
%R  " ^ 
^MA  " 


Then 


timnel  radius  = A - A-r-  x CLIN  = 2.34  x 10“-^  meter 

*2 

radial  mesh  = .635  x 10“"^  meter 

axial  mesh  x CLIN/N.,.  = .9525  x 10“^  meter. 


Coil  parameters  converted  to  MKS  units: 


Radius 

Axial  Pos. 

Strength 

#1 

3.7592x10"^ 

-3.81x10"^ 

-4.0181x10"^ 

#2 

3.7592x10"^ 

3.81x10"^ 

4.3129x10"^ 

#3 

3.7592x10"^ 

11.43x10”^ 

-4.3129x10"^ 

#4 

3.7592.10"^ 

19.05x10"^ 

4.0183x10"^ 

We  will  hand- calculate  the  matrix  entry  for  the  poinu  one 
mesh  \mit  off  the  axis  (R  = .635  x 10~5)  at  the  center  of 
the  first  gap  (2  = 2^1=  3.81  x 10~^):  the  dominant  contribu- 
tion is  from  coil  #2,  for  which 

c = 4x. 635x10"^  X 3.7592x10"^/  |(3. 7592x10"^ 

+ .635x10"^)^  + 0^1]  = .703 


TvkTs 


Figure  10:  Geometry  for  Magnetic  Vector  Potential  at 

(R,Z)  Due  to  Current  I in  Loop  Radius  R^ 
Centered  at  (0,Z^) 


To  obtain  K and  E from  tables,  which  are  usually  given  in 
terms  of  the  modulus  expressed  as  an  angle,  we  take  arcsin 
.703  = 44.7°  as  our  enti^.  In  the  program,  subroutine  ELIVA 
uses  the  parameter  m = as  input.  By  either  method,  we  ob- 
tain K(c)  = 1.8493,  E(c)  = 1.3535,  and 

M2  = 4.3129x10“^  (.635x10"^  x 3-7592x10"^)^^^ 

= 7.361x10"® 

For  coil  #3,  Z = 3.81x10"^,  Z-|  = 11.43x10"^,  giving  c = .3513, 

K = 1.6229,  E = 1.5212,  M3  = -.591x10“8.  Similarly  the  con- 
tribution from  coil  #1  is  M-)  = -.550x10~8.  pcr  coil  #4,  c = 

.1948.  Since  this  is  less  than  .2,  we  use  (3-9)  to  evaluate 
the  square  bracket  in  (3-5): 

M^  = 4.0183x10"^  (.635x10"^  x 3.7592x10"^)^^^ 

X (1  + (3/4)(.1948)^  + (.1948)^j  = .093x10"® 

Thus  the  resultant  M is 

(7. 361  - .591  - .550  + .093)  X 10"®  = 6.31x10"® 


Subroutine  MAGVA  carries  out  this  computation  for  each  of  the 
mesh  points  defined  by  (3-7)  and  (3-8).  Ihe  run  for  the  stan- 
dard case,  using  the  coarse  4 x 16  matrix  size  to  fit  on  the 
time-sharing  system,  is  shown  in  Fij^reH.  The  complete  vector 
potential  matrix,  multiplied  by  10'^  and  rounded  to  integer 
values  for  clarity,  is  given  in  reverse  row  order. 

The  peak  value  623x10  on  the  row  next  to  the  axis  (row  of 
zeroes)  corresponds  to  the  6.31x10“^  we  have  just  hand-calcu- 
lated. The  discrepancy  is  due  to  the  fact  that  hand  calcula- 
tion with  linear  interpolation  of  tables  is  barely  adequate 
for  this  problem;  it  does  serve  as  a useful  check  that  no  gross 
errors  have  been  made.  The  machine  calculated  values  are  cor- 
rect to  4 or  more  digits. 

The  axial  and  radial  fields,  in  gauss,  are  tab\ilated  below  the 
matrix. 
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A further  check  on  the  correctness  of  a matrix  of  this  kind 
(which  should  be  applied  as  a check  if  any  different  method 
of  computing  the  matrix  is  adopted) , is  obtained  by  examining 
any  entry  on  the  row  adjacent  to  the  axis.  The  field  at  the 
corresponding  point  on  the  axis  is  then 


B 


gauss 


2x10^° 


x (matrix  value  adjacent  to  axis) 

p 

(radial  mesh  size  in  mm) 


(3-10) 


Thus  for  the  point  calcialated,  using  the  machine  value  for 
the  off-axis  potential t 


B ^ 2x10^^  X 623x10'^^ 

gauss  .635^ 


3090 


This  is  in  agreement  with  the  value  shown  for  this  point  in 
the  axial  field  tabulation  in  Figure  8,  but  is  1.4%  higher 
than  the  expected  3046.7  gauss.  This  discrepancy  is  not  due 
to  inaccuracy  of  computation,  but  is  simply  a result  of  the 
coarse  4 x 16  mesh  used  in  this  excinple  — it  is  a discreti- 
zation error.  If  the  same  data  is  run  with  20  meshes  radially 
instead  of  4,  the  computed  peak  axial  field  is  3049  gauss,  an 
error  of  less  than  0.1%.  This  point  should  be  remembered  if 
(3-10)  is  used  to  check  any  alternative  method  of  computing 
MAGMAT. 

Since  actual  fields  in  tubes  are  seldom  known  to  better  than 
5%  accuracy,  this  also  indicates  that  the  discretization  error 
is  not  serious,  in  practical  terms,  even  for  the  coarse  4 rad- 
ial mesh  case. 


4.0  THE  RF  VECTOR  POTENTIAL  MATRIX 


This  matrix  provides  the  working  description  of  the  rf  fields 
in  the  tunnel  and  gap  region.  Its  computation  thus  depends  on 
the  boundary  conditions  that  are  assumed. 


4. 1 Basic  Field  Model 


Following  the  lead  of  Kosmahl  and  Branch  we  adopt  a field 
model  in  which  the  gap  field  increases  in  intensity  nerr  the 
noses,  but  not  to  the  extreme  values  associated  with  a sharp 
edge.  Kosmahl  and  Branch  take  the  axial  gap  field  at  the  tun- 
nel radius  to  vary  as  cosh(mz)  where  z is  measured  from  the 
gap  center,  and  m is  an  arbitrary  parameter  (^0)  which  in 
effect  describes  the  'sharpness’  of  the  noses.  K.  and  3.  give 
experimental  data  confirming  that  the  model  is  a good  one  for 
a typical  nose  radius,  giving  a field  concentration  of  about 
2.5:1  at  the  nose  compared  to  the  gap'  center.  This  corresponds 
to  ml  = arc  cosh  2.5  = 1.57,  where  I is  the  half  gap  length. 

It  should  be  noted  that  the  model. does  have  a logical  incon- 
sistency, in  that  the  finite  field  concentration  corresponds 
to  a rounded  nose,  but  a rounded  nose  does  not  correspond  ex- 
actly to  a boundary  condition  of  = cosh(mz)  up  to  i and  zero 
beyond,  since  the  tunnel  wall  curves  up  before  it  reaches  z = t. 
Thus  the  model  will  break  down  if  one  tries  to  examine  the 
fields  near  the  nose  on  a scale  comparable  with  the  implied 
nose  radius,  or  to  specify  an  excessively  large  field  concen- 
tration factor  — about  4 is  the  limit  that  should  be  used, 
and  1.5  to  3 is  a more  reasonable  range.  This  is  in  general 
agreement  with  the  conclusions  of  K.  and  B.  The  case  m=0, 
concentration  factor  1 , corresponds  to  the  uniform  gap  f.ield 
first  analyzed  by  Wang  [lOl.  The  program  will  call  for  the 
field  concentration  factor  as  input,  and  will  calculate  m from 
it  (knowing  t),  because  it  is  easier  for  the  user  to  think  in 
terms  of  concentration  factor. 


[9J 


H.G.  Kosmahl  and  G.M.  Branch,  'Generalized  Representation 
of  Electric  Fields  in  Interaction  Gaps  of  Klystrons  and 
Traveling  Wave  Tubes' , IEEE  Transactions  on  Electron 
Devices,  ED-20  #7,  July  1973,  pp.  621-629. 


C.C.  Wang,  'Electromagnetic  Field  Inside  a Cylinder  with 
a Gap',  Journal  of  Applied  Physics,  1_6,  June  1945, 
pp.  351-366. 
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The  mid-gap  field  strength  (Eq)  and  the  total,  rf  voltage 
across  the  gap  (V^f)  are  related  by 


^0  2 sinh  (mt) 

If  m=0,  this  reduces  to  as  expected.  We  are  going 

to  compute  the  matrix  for  Vj.f=1  volt  peak. 


• 4.2  Vector  Potential  Expressions 

It  will  be  remembered  from  the  introduction  tliat  the  quantity 
Vec  to  be . stored  in  matrix  is  ' radius  x vector  potential ' , 
just  as  in  the  magnetic  field  case.  Mioltiplying  the  four  ex- 
pressions given  by  Kosmahl  and  Branch  (for  the  axial  and  radial 
fields  in  the  gap  and  in  the  tunnel)  by  r,  and  integrating  with 
respect  to  r and  z,  we  obtain  the  following  two  expressions  for 
Vec.: 

For  C ^\z\^  1: 

rJ..^(r/k^i-m^) 

^+m^  (a'/k^+m^) 


= K 

8C  0 


cosh(mz) 


- a 


p^-ma  Pj^+ma 


cosh 


Pn^ 


(4-2) 


and  for  z ^ i: 


V = E 
ec  c 


rJ,(k  r) 


V In,  ■ 

4,  Pn^l^V*  1 

n_i 


sinh(Pj^t/a+mt)  sinh(p^t/a-mt) 


Pjj+ma 


p -ma 


>e 


-Pn|z|/a 


('•-5) 
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wherein 


Eq  = . field  (v/m)  at  z=0,  r=a,  defiixed  by  (4-1) 
a = tunnel  radius,  meters 
I = half  gap  length t meters 
(JO  = angular  frequency 
k = w/c 

Xjj  = nth  root  of 

k = X /a 
n n 

The  geometry  is  illustrated  in  Figure. 

It  may  be  verified  by  differentiating  (4-2)  and  (4-3)  with 
respect  to  r and  z,  using  the  known  relation  [11] 


.J^(kx))  = kxJ^(kx), 


that  the  results  agree  v;ith  K and  B' s expressions  for  the 

field  components,  including  the  factor  1/r  which  we  need.  i 

Thus  with  V defined  by  (4-2)  and  (4-3),  we  have  for  the 

fields 


E„  = 


1 

ec 


E„ 


1 


“5F 


(4-4) 

(4-5) 


The  negative  sign  in  (4-5)  will  be  taken  care  of  later. 

The  units  of  Vg^  are  meter- volts,  and  and  E^,  will  then 
be  in  volts  per  meter. 

The  matrix  VgQ  is  symmetrical  about  z=0,  so  the  left-hand 
half  (z  < 0)  can  be  filled  in  once  the  z ^ 0 values  have 
been  calculated  (assuming  that  we  were  not  so  stupid  as  to 
adopt  an  unsymmetrical  grid) . 


[11J  N.W.  McLachlan  'Bessel  Functions  for  Engineers',  Oxford 
Univ.  Press,  p.  158,  equations  22  and  24. 
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4.3  Ccroputation  of  the  rf  Vector  Potential  Matrix 


We  shall  compute  one  matrix  Vec  for  unit  rf  voltage  across 
the  gap;  this  describes  the  shape  of  the  field  for  any  cavity. 
When  a super  electron  is  in  cavity  n with  peak  rf  voltage  V„ 
at  phase  0^*  'conceptually  lay  Vgc  over  that  cavity,  obtain 
and  at  the  position  of  the  super  electron  by  (4-4)  and 
(^5),  and  multiply  by  Vncos(t»t+0n) . Thus  one  matrix  serves 
for  all  cavities. 

Ihe  immediate  problem  is  thus  to  devise  an  efficient  scheme 
for  evaluating  (4-2)  or  (4-3)  at  each  node  of  a grid  of  the 
desired  fineness,  as  shown  in  Figure  12..  The  constajjt  factor 
Eq,  given  by  (4-1)- with  can  be  omitted  for  the  present. 

We  find  that  in  some  regions,  50  or  more  terms  are  needed  to 
get  5 figure  accuracy,  while  in  other  regions  less  than  ^0  are 
adequate.  Thus  the  computation  scheme  does  warrant  some  care- 
ful thought,  because  50  terms  of  (4-2)  or  (4-3)  is  obviously 
not  a trivial  calculation.  But  we  can  afford  to  be  fairly 
generous  in  the  number  of  terms,  because  we  only  calculate  one 
matrix  one  time;  compared  with  direct  calculation  of  the  fields 
at  each  superelectron  position  at  each  time  step,  the  vector 
potential  approach  will  start  to  return  dividends  in  time  saved 
after  only  about  2 time  steps  of  the  main  calculation. 

The  complexity  of  (4-2)  and  (4-3)  also  makes  it  unlikely  that 
we  could  establish  analytically  the  number  of  terms  required 
for  a given  accuracy  in  a given  region  of  r and  z.  So  we  pro- 
ceed heuristically  by  summing  the  series  directly  to  a consid- 
erable number  of  terms  at  some  representative  points  both  in 
the  gap  and  in  the  tunnel,  printing  out  the  partial  sum  after 
each  term.  This  is  dene  for  a 'typica].*  gap  to  diameter  ratio 
and  a 'typical'  diameter  to  wavelength  ratio.  Inspection  of 
the  output  establishes  for_each  point  the  number  of  terms  needed 
to  get  within  1 part  in  10^  of  the  final  potential.  The  results 
of  such  a computation  are  shown  in  Figure  1Z.  We  see  that  in 
the  middle  2/3  of  the  gap,  and  in  the  tunnel  beyond  about  4/3 
of  the  half  gap  length,  8 terms  or  less  are  sufficient.  As  we 
approach  the  plane  z=t  (the  gap  edge)  from  either  side,  the 
number  of  terms  rises  first  to  about  16,  then  to  about  30,  and 
at  z=l  and  r=a  even  60  are  inadequate.  Taking  a generous  mar- 
gin to  allow  for  different  frequencies  and  l/a  ratios,  we  ar- 
range to  sum  12  terms  for  z/t  < .7  or  z/t  > 1.3,  and  20  terms 
for  .7  ^ z/t  ^ .93  or  1.07  ^ z/l  ^ 1.3.  The  region  .93  < z/l 
< 1.07  is  clearly  one  of  slow  convergence  which  will  require 
different  treatment. 
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The  summations  to  12  or  20  terms  are  carried  out  by  direct 
evaluation  and  summation  of  the  terms  as  written.  The  sub- 
routine BESVA  is  used  for  evaluation  of  the  J-](kjjr)  terms; 
it  is  more  than  twice  as  fast  as  the  Library  routine  BESJ> 
and  changes  over  automatically  to  the  trigonometrical  ex- 
pansions when  k„r  > 12.  The  J-jCkna)  terms  do  not  have  to 
be  calculated;  oy  definition,  kna=Xj^  which  is  a root  of  Jq, 
and  the  corresponding  J>]  values  are  tabulated;  the  table  up 
to  n=20  is  incorporated  in  the  program,  with  the  correspond- 
ing table  for 


4.4  Special  Methods  near  the  Gap  Edge 

Figure  shows  the  convergence  for  two  po.  nts  near  the  plane 
z=t;  the  upper  curve  at  z=.9t,  r=.9a  shows  that  20  terms  were 
adequate  there,  but  the  lo\/er  curve  at  z=.95'^  shows  that  about 
30  terms  should  be  used  here.  As  a function  of  r,  we  find 
that  the  convergence  is  more  rapid  for  middle  values  of  r, 
but  is  slowest  for  r/a  ^ .1  or  ^ .9.  The  region  near  the 
axis  is  less  important,  because  only  a small  part  of  the  beam 
travels  there,  so  we  concentrate  on  the  high  values  of  r/a, 
but  less  than  unity.  Figure  J3f<»Xagain  for  z=.9t  and  .95t, 
but  at  r=a,  shows  that  tha  oscillatory  nature  of  the  conver- 
gence has  now  disappeared.  But  it  is  still  true  that  20 
terms  are  adequate  at  z=.9t,  while  30  or  so  are  needed  at 
.95^. 


Beyond  n=10,  we  can  start  to  make  simplifying  approximations 
before  continuing  the  summation,  because  ^ and  p^  are  now 
both  greater  than  60;  for  example,  coshCp^z/a)  in  (4-2)  can 
be  replaced  by  .5  exp(pj^z/a),  and  this  can  then  be  combined 
with  the  exp ( -PnVa ) term  to  give  exp(-pjj(t-z)/a) . With 
Pn  > 60,  this  is  going  to  give  rapid  convergence  except  when 
z « t — now  v/e  see  why  the  region  is  the  most  difficult 
part. 


4.5  Approximations  Valid  for  Large  n 

Specifically,  the  approximations  we  adopt  for  n > 20  are 


= TT(n-.25) 

(4-6) 

= TT(n-.25)/a 

(4-7) 
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Tr 


(4-8) 


cosh(py^z/a)  = .5  exp(Pj^z/a)  (4-9) 

‘ - i]  (^10) 

except  when  <12 

sinh(pj^t/a+mt)  = .5  exp(pj^t/a+mt)  (4-11) 


We  note  that  ^ within  .04?6,  but  since  the  square  root 

in  the  exact  definition  of  Pj,  is  a very  fast  operation,  we 
need  not  take  this  approximation. 

The  approximations  (4-6)  and  (4-8)  are  better  than  .0^%  even 
at  n=12,  so  we  are  quite  safe  in  adopting  them  for  n > 20. 

The  limitation  k^r  < 12  will  only  affect  the  innermost  rows, 
if  any;  at  n=21,  kj^=65.2/a,  so  the  limitation  is  equivalent 
to  r/a  < .185,  which  applies  to  no  row  of  the  matrix  if 
Ncr  ^ 5,  and  only  one  rov/  for  up  to  10,  which  covers  most 
cases;  as  stated  ir.  the  introduction,  we  do  not  expect  Nqr  ever 
to  exceed  20.  For  n > 21,  the  limitation  becomes  progressively 
less  significant.  Applying  these  approximations  and  simplify- 
ing, we  find  that  the  general  terms  in  the  summations  for  both 
(4-2)  and  (4-3)  can  be  expressed  as  . 


Var 

— sin 
^n 


^n  TT 


exp(-p^|t-z|/a)  F(a,m,t,n) 


(4-12) 


where  F(a,m,t,n)  = - i ^ 


p -ma  p„+„„ 
*^n  -^n  ma 


for  \z\  ^ I 


(4-13) 


mt  -mt 

and  F(a,m,i,n)  = for  |z|a  I 

n -^n 


(4-14) 


for  n > 20. 
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Using  these  approximations,  we  continue  the  summation  out  to 
n=40,  for  .93  < z/t  < .99  or  1.01  < z/t  <1.07.  We  shall  make 
a further  adjustment  of  n after  the  next  section. 


4.6 


The  Case  z=t 


At  z=t,  we  find  that  40,  or  even  60  terms  are  not  adequate, 
and  we  look  for  a more  sophisticated  approach.  This  case  is 
not  an  improbable  one,  because  of  the  tendency  to  choose 
'round  numbers'  for  setting  up  cases.  For  example,  if  the 
cavity  period  is  0.3"  and  the  gap  C.1",  any  choice  cf  the  num- 
ber of  axial  meshes  Nqa  that  is  a multiple  of  6 will  place  a 
grid  line  exactly  on  ^e  plane  z=t.  Remember  that  we  are  not 
required  to  compute  for  arbitrary  values  of  z,  but  only  for 
those  corresponding  to  a line  of  the  chosen  grid.  Further, 
a distinction  betv/een  2.1 1 - .99  and  z/t  = 1.00  is  not  very 
meaningfifL  in  terms  of  typical  TWT  dimensions  and  achievable 
tolerances,  so  we  shall  take  z=t  if  the  nominal  z is  between 
.991  and  1 .01 t. 

On  this  plane  we  can  use  either  (4-2)  or  (4-3),  and  we  should, 
of  course,  get  the  same  result.  Figure  13(c)  shov/s  the  conver- 
gence at  2=1,  r=.9a,  calculated  both  ways.  Clearly  both  curves 
are  converging  to  a value  of  about  1.746x10“°,  but  have  not  con- 
verged within  an  acceptable  tolerance  even  at  50  terms;  by 
chance,  this  happens  to  be  a particular  number  of  terms  at 
which  they  both  cross  over  the  asymptotic  value,  as  are  41 
or  31  terms.  The  periodicity  of  the  curves,  and  hence  the 
specific  favorable  values  of  n,  depends  on  the  ratio  r/a.  The 
oscillatory  component  comes  from  the  sin  term  in  (4-12),  and 
clearly  we  can  obtain  satisfactory  accuracy  without  an  exces- 
sive number  of  terms  if  we  stop  at  one  of  the  crossovers. 

4.6.1  Diophantine  Approximation  for  r<9,  z=t 

Since  the  curve  is  effectively  the  integral  of  (4-12),  it  is 
the  cosine  function  that  should  be  zero;  thus  we  should  choose 
n to  make  j an  integer  (or  almost  so)  in 


1 

I 


I 

a 


Pn^  TT 


(o+i)^ 


Using  the  approximation  p^  = V*  " <r(n-.25),  we  have 


n = .25+(a/r)( 0+.75) 
or  0 = (n-.25)(r/a)“.75 
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(4-15) 

(4-16) 


0 


10 


20  30 


60 


Figure  13(c) 


The  summation  will  then  be  terminated  close  to  one  of  the 
asymptote  crossovers  rather  than  at  one  of  the  peaks.  Since 
n must  be  an  integer,  and  j should  be  as  close  to  an  integer 
as  possible,  this  leads  us  more  or  less  into  the  realm  of 
Diophantine  [12 j equations  (algebraic  equations  restricted 
to  integer  solutions),  which  are  notorious  for  not  having 
any  general  methods  of  solution-  Knowing  this  to  be  the 
case,  we  shall  evaluate  (4-12)  out  to  some  fixed  n,  say  40, 
chosen  to  get  dov/n  to  the  accuracy  region,  and  then  on 
to  the  next  integer  n satisfying  (4-15);  we  shall  obtain  this 
by  solving  (4-16)  to  find  a j that  is  acceptably  close  to 
being  an  integer  — in  general  there  v;ill  not  be  a strict 
Diophantine  solution  except  for  some  particular  values  of  r/a. 

Subroutine  EANTUS  identifies  the  next  crossover  for  any 
given  starting  value  of  n and  r/a,  excep+.  r/a=1;  if  it  hap- 
pens that  n is  itself  a. crossover  point,  FANTUS  fails  to 
recognize  this  and  goes  out  to  the  next  crossover,  but  other- 
wise it  finds  the  first  available  one.  For  Nqr  ^ 20,  r/a 
will  never  be  less  than  .05  or  greater  than  .95,  and  vre  find 
that  the  majamum  number  of  extra  terms  called  for  by  FANTUS 
is  20;  for  the  more  likely  value  of  8 for  a maximum  of  9 

extra  terms  is  needed,  and  in  the  middle  range  of  r/a  values 
it  is  down  to  5 or  less.  For  small  r/a,  and  z/t  close  to  1, 
the  convergence  pattern  is  of  the  'beating  wave'  form  shown 
in  Figure  14.  Starting  from  an  arbitrary  point  such  as  A, 
FANTUS  correctly  identifies  the  next  envelope  crossover  at 
3,  and  is  not  deceived  by  the  intermediate  point-to-point 
crossovers . 

Since  the  oscillatory  term  in  (4-12)  is  independent  of  z/4, 
this  theory  is  equally  valid  for  optimizing  the  number  of 
terms  near  z=t  as  well  as  on  it,  so  we  apply  it  over  the  whole 
range  covered  by  Section  4.5,  even  though  it  may  not  be  strict- 
ly necessary  there.  Thus  for  .95  ^ z/t  < -99  and  1.01  < z/t 
^ 1.07  for  all  r/a,  and  for  .99  ^ a/t  ^ 1.01  for  r/a  ^ .95,  we 
extend  the  summation  of  (4-12)  beyond  40  terms  out  to  an  opti- 
mum number  between  40  and  60  determined  by  FANTUS  (always  60 
for  r=a).  The  values  determined  by  FANTUS  for  Nrp^  = 20  are 
tabulated  in  Figure  17  as  a function  of  the  radial  mesh  number 
I. 


[^12]  Diophantus,  ' Arithmetica' , Univ.  of  Alexandria,  Egypt, 
ca.  320,  trans.  S.  Stevin,  pub.  Elsevier,  Leyden,  1634. 
(Newton  collection,  Bender  Library,  Stanford  University) 
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4.6.2  2eta  Function  Approximation  foi‘  r=a,  z=-t 


At  r=a,  the  Diophantine  equation  (4-16)  has  no  solution,  be- 
cause the  oscillatory  term  in  the  expansion  has  disappeared. 
Figure  V?  shows  that  the  convergence,  calculated  from  either 
side,  is  now  monotonic,  but  so  slov;  that  even  60  terms  are 
quite  inadequate.  But  at  n=60,  a further  simplification  can 
be  made:  ma  is  now  negligible  compared  to  Pj^  and  the  residue 
of  the  summation  reduces  to 


2 

-a  CO 


sh(mt)  ^ 
61 


cosh(mt) 


0> 


I 


1 

(n-.25)^ 


(4-17) 


This  is  a Riemann  Zeta  function  in  the  generalized  form 
introduced  by  Hurwitz.  Tables  of  the  generalized  functions 
are  not  readily  available,  but  the  sum  can  be  reduced  to 
known  forms  as  follows: 


61 


~ - 60 

Z XTT-;25)2  = Z (H-^25)^  ■ Z 
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= ^ Xn-^2572  “ 2.5252825  by  direct  summation 
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Now  the  first  summation  is  Riemann's  c(^)  = 1.64493^7  [13] 
and  the  second  is  Catalan's  constant  X ==..91596559  [I3f  p.  80?1 
Hence  " 

00 

y 7 - .0165972*  (4-18) 

^ (n-.25)2 

Hence  the  residue  of  the  summation  is 

a2 

RgQ  = - .0165972  “2  cosh  ml 

V 

= - 1.682  E-3  s?  cosh  mt  (4-19) 


Note  that  the  nxjmerical  coefficient  1.682E-3  is  specific  to 
stopping  the  term-by-term  summation  at  60  terms.  We  adopt 
this  formulation  if  z is  within  of  equality  with  i and 
r=a.  Outside  this  range,  the  terms  equivalent  to  e”Pn 
ensure  convergence  within  60  terms  for  any  reasonable  t/a. 

The  result  of  summing  from  the  gap  side  to  60  terms  and  then 
adding  the  (negative)  residue  R50  is  shown  by  the  x in  Figure 
15.  Clearly  it  has,  in  this  instance  at  least,  hit  the  aver- 
age of  the  two  curves  very  closely,  while  only  requiring  one 
series  to  be  summed. 

Figure  16  shows  the  variation  of  Vec  at  r=a  going  through  the 
z=t  region,  indicating  that  the  various  methods  used  do  join 
up  smoothly. 

To  summarize,  the  computation  strategy  is; 

For  z/t  < .7f  sum  (4-2)  to  12  terms. 

For  .7  ^ z/t  < .93,  sum  (4-2)  to  20  terms. 

For  .93  ^ z/t  < .99,  for  all  r/a,  and  for  .99  ^ z/t  ^ 1.01  for 
r/a  ^ .95,  sum  (4-2)  to  20  terms,  then  (4-12)  to  (40  to  60) 
terms  as  determined  by  FANTUS. 

For  .99  ^ z/t  ^ 1.01  and  r=a,  sum  (4-2)  to  20  terms,  then  (4-12) 
to  60  terms,  then  add  (4-19). 


mj  M.  Abramowitz  and  I.  Stegun,  'Handbook  of  Mathematical 
Functions',  N.B.S.  Washington,  D.C.  1964  or  Dover  Publi- 
cations, New  York,  1965,  page  811. 

*can  be  obtained  more  briefly  by  Gumowski's  method,  J.A.P.  August 
1953,  p.1068  (with  correction  noted  on  p.  1330).  This  gives 
.0165971.  We  did  not  find  this  reference  until  after  this  report 
was  first  issued. 
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Figure  16 


For  1.01  < zA  ^ 1.07,  sum  (4-3)  to  20  terms,  then  (A-12)  to 
(40  to  60)  terms. 

For  1.07  < zA  ^ 1.3,  sum  (4-3)  to  20  terms. 

For  z/t  >1.3,  sum  (4-3)  to  12  terms. 


4 . 7 Subroutine  LALAVA 


We  have  now  established  procedures  for  evaluating  (4-2)  and 
(4-3)  with  the  necessary  accuracy  in  all  regions  of  interest. 
Subroutine  LALAVA  carries  out  these  evaluations  at  all  nodes 
of  the  chosen  grid  for  z > 0,  and  copies  them  to  the  corres- 
ponding z < 0 nodes.  The  values  on  the  axis  are  all  zero, 
without  computation-,  and  the  values  at  R = -1  are  equal  to 
those  at  R = 1 . The  Nqr+I  row  is  equal  to  the  Nqj^-I  row  for 


|z  I > t,  thus  forcing  the  tangential  component  to  be  zero 


at  the  wall.  In  the  gap  region,  the  Njqj^+I  row  is  extrapolated 
from  the  and  Nq^^-I  rows  to  maintain  the  required  axial 


fields. 


As  in  the  case  of  the  magnetic  vector  potential  matrix,  the 
options  are 


i) 

ii) 


compute  the  matrix,  use  it  and  discard  it. 


compute  matrix  and  store  on  file  REMAT,  as  well  as 
using  it  for  the  current  run. 


iii)  read  in  the  matrix  from  RFMAT. 


When  the  matrix  is  stored,  it  is  preceded  by  the  nominal  dimen- 
sions Ncr  and  Nq^,  as  a safeguard  against  reading  in  the  wrong 
matrix. 


LALAVA  will  also  print  out  the  matrix  if  desired,  and  will  also 
compute  the  field  components  (4-4)  and  (4-5)  at  each  point  and 
print  out  tables  of  Ej,  and  E^.  Figures  16,  17  and  18  are  ex- 
amples of  these  printouts,  only  half  the  region  being  shown  in 


each  case, 
symmetric. 


Vec  and  E^, 


are  symmetric  about  z=0  and  E^  is  anti- 


As  a final  check,  the  fields  were  also  computed  directly  from 


Kosmahl  and  Branch's  expressions,  and  the  result  is  shown  in 


Figure  19.  The  comparison  of  Figure  19  with  Figure  18  is  not 
a precision  one,  because  we  did  not  go  through  -the  whole  rou- 
tine of  finding  appropriate  large  n simplifications  for  each 
region,  but  simply  summed  20  terms  at  every  node.  The  agree- 
ment is  excellent  in  the  mid-gap  region,  and  inside  the  tunnel, 
but  degrades  near  the  gap  edge,  as  woiild  be  expected  with  only 
20  terms  taken.  The  field  expressions  have  even  slower  conver- 
gence than  the  potential  expressions,-  so  that  still  more  terms 
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Figure  17 
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Figure  18 


wolild  have  to  be  taken  to  get  a precision  comparison;  the 
comparison  shown  is  good  enough  to  demonstrate  that  no  mis- 
takes in  scaling  have  been  made;  it  also  shows  that  the  small 
negative  values  of  for  large  r and  small  but  non-zero  z, 
which  are  not  realistic,  are  a basic  defect  of  the  model,  not 
an  effect  of  using  the  potential  method.  It  is  a side  effect 
of  the  inconsistency  of  the  model  at  the  gap  edge,  discussed 
earlier. 

Despite  this  deficiency,  which  is  numerically  not  very  large, 
it  is  the  opinion  of  this  writer  that  the  K and  B model  is 
the  best  one  to  which  we  know  an  analytic  solution.  The  only 
model  v/hich  can  in  principle  deal  correctly  with  noses  of  fin- 
ite radius  is  the  relaxation  method  on  a siiitably  fine  mesh. 
This  is  the  approach  used  in  the  Los  Ai.amos  program  LALA;  the 
reason  for  not  regularly  using  LALA  is  a matter  of  size  and 
time : LALA  has  about  60  pages  of  source  statements  if  close- 
packed  (actually  94  pages  as  normally  printed),  and  typical 
solution  times  are  300  to  600  seconds  of  CPU  time.  The  ana- 
lytic subroutine  LALAVA  developed  from  the  foregoing  analysis 
occupies  5 pages  of  source  statements,  and  has  tjrpical  solu- 
tion times  of  20  to  40  seconds.  However,  the  option  for  read- 
ing in  a previous  LALAVA  matrix  will  be  written  so  that  a ma- 
trix generated  by  LALA  would  also  be  accepted;  there  will  be 
problems  of  adjusting  the  scale  factor,  since  LALAVA  normal- 
izes to  unit  rf  voltage  across  the  gap,  which  is  the  important 
quantity  for  TWT  work,  while  LALA  normalizes  with  respect  to 
energy  chonge  along  the  axis,  which  is  the  important  quantity 
for  accrl  ''■ators,  for  which  LALA  was  originally  written. 

In  Sectic)  2.1  't  was  asserted  that  the  fields  at  the  mid- 
planes of  che  tunnels  woui d be  25  to  30  dB  below  the  gap 
fields,  so  that  beyond  these  planes  the  fields  could  safely 
be  neglected.  Figure  18  shows  that  the  field  on  the  axis  at 
the  tunnel  mid-plane  (7.0  volts/meter)  is  below  the  mid- gap 
axial  field  (195.4  v/m)  by  28.9  dB,  and  is  falling  by  a fur- 
ther 2.5  dB  per  mesh  point.  Thus  the  assertion  is  well  sup- 
ported in  this  numerical  case,  which  is  a quite  typical  one. 


5.0  SPACE  CHARGE  FORCES 


As  in  the  case  of  rf  and  magnetic  fields,  the  space  charge 
forces  are  to  be  derived  from  the  gradients  of  a matrix  of 
potentials.  This  matrix  differs  from  the  others  in  that  it 
is  moving  with  the  beam,  and  that  it  has  to  be  recalciilated 
completely  after  every  time  step  of  the  trajectory  calcula- 
tion, since  the  distribution  of  space  charge  changes  at  each 
step.  For  this  reason,  the  most  extreme  efforts  must  be  made 
to  obtain  a fast  and  efficient  algorithm  for  this  matrix. 

For  a rectangular  geometry,  the  fastest  known  numerical  solu- 
tion of  Poisson's  equation  is  the  Hockney- Buneman  FACR 
(Fourier  Analysis  Cyclic  Reduction)  method  What  fol- 

lows is  primarily  an  extension  of  this  method  to  cylindrical 
coordinates. 

The  essence  of  the  FACR  method  is; 

i)  assignment  of  the  continuous  distribution  of  charge 
into  discrete  charges  at  the  nodes  of  a mesh. 

ii)  Fourier  analysis  of  ttie  charge  d.‘ stribution  in  one 
direction,  along  each  row  of  the  mesh. 

iii)  combination  of  the  analyzed  rov/s  in  sets  of  3,  using 

certain  trigonometric  identities  to  eliminate  alternate 
rows,  so  that  the  number  of  rows  left  is  reduced  by  a 
factor  2. 

iv)  Repeating  this  cycle  until  it  is  reduced  to  a relation 
between  the  center  row  and  the  boundary  rows,  on  v/hich 
the  potentials  are  determined. 

v)  Reversing  the  process  to  fill  in  the  potentials  on 
intervening  rows  in  the  reverse  order. 

The  difficulties  encountered  in  applying  this  technique  in 
cylindrical  coordinates  are  two; 

i)  The  expressions  for  cyclic  reduction  are  deper.dent  on  r, 
whereas  the  corresponding  rectangular  geometry  expres- 
sions were  independent  of  y. 


[14J  R.W.  Hockney,  'The  potential  Calculation  and  some 

Applications',  in  'Methods  in  Computational  Physics', 
Ed.  B.  Alder  et  al.  Academic  Press,  New  York  1970. 


ii)  One  of  the  boundaries  is  the  axis,  on  which  the  poten- 
tials are  not  initially  known.  Instead  we  have  the 
condition  that  equipotentials  must  intersect  the  axis 
at  right  angles. 

In  this  section  -we  shall  develop  the  explicit  recursion  rela- 
tions for  the  cylindrical  case,  and  show  how  to  handle  the 
axis  by  developing  a second  recursion  to  be  solved  simultan- 
eously with  the  first. 


5.1  Analysis  of  the  Axisymmetric  Space  Charge 
We  require  the  solution  of  Poisson's  equation 
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r Ir 
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(5-1) 


at  the  nodes  of  the  grid  shown  in  Figure  20.  The  potential  P 
is  in  this  case  the  ordinary  scalar  potential.  p(r,z)  is  the 
charge  density  at  r,z  in  coulombs/neter^;  it  will  be  determined 
from  the  superelectrons  in  the  vicinity  of  (r,z)  by  formulas  to 
be  given  later. 

Longitudinally  the  g^id  extends  over  one  beam  wavelength  Xg 
plus  one  mesh,  from  J=0  to  J=Nc^+1 . Nsa  is  one  of  the  numbers 
6,  12,  24  or  possibly  48,  for  which  very  fast  Fourier  transform 
routines  exist.  The  mesh  length  is  thus 


hsa  = (5-2) 

In  the  radial  direction,  Nsr  must  be  a power  of  2 to  allow  the 
Cyclic  Reduction  to  come  do\m  to  a single  row  half  way  between 
the  axis  and  the  wall.  Usually  Nsr  will  be  4 or  8,  possibly  2 
or  16.  We  define  a shape  factor  f such  tliat  the  radial  mesh 
size  is  fh  and 


h„^  = (fhor>  = s/Ncp 
sr  sa  SR 


(5-3) 


The  charge  Q(i,j)  to  be  associated  with  node  (i,j)  of  the  grid 


Q(i»0)  =— T-^p(r,z) 


(5-4) 
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The  node  is  located  at  r=ifh,  z=Jh. 


The  differentials  in  (5-1)  are  now  expressed  in  finite 
difference  form  [^15]]: 


Sr 


= (r(i-i,j)  - 2p(i,j)  + p(i+i,o)| 


(5-5)  L 


I?  = m |p(i+^.0)  - P(i-1,0)j 


(5-6) 


^ (p(i-i,j)  - 2P(i,o)  + P(i+1,o)] 

Sz  h 1 : J 


(5-7) 


Substituting  (5-5)  throii.gh  (5-7)  in  (5-1) 


(1+1/2i)  P(i+1,o)  + (1-l/2i)  P(i-1,o)  + f^|p(i,o-1)  + P(i,o+1)j 


- 2(1+f^)  P(i,o)  = Q(i,o) 


(5-8) 


for  i = 1 to  Nsr-1. 
On  the  axis; 


2P(1,j)  + f2|p(0,j-1)  + P(0,o+1)|-  2(l+f^)  P(0,j)  = Q(0,j)  (5-9) 


We  have  now  split  Poisson's  equation  into  Nqr  separate 
equations,  which  differ  because  of  the  1-l/2i  and  1+1 /2i 
factors. 


L15J  F.S.  Shaw’  'Relaxation  Methods',  Dover  Publications, 
New  York,  1953. 
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5.2  Fourier  Analysis  Step 


We  now  express  the  unknown  potential  P and  the  known  or  given 
charge  distribution  Q on  each  row  of  the  grid,  as  finite 
Fourier  sums.  These  sums  extend  to  the  mth  harmonic,  where 

m = 5 Ngj^  (5-10) 


m 

m-1 

P(i,j) 

= U(i,k)  cos  ^ 0 + 

Y.  ^ ^ 

(5-11) 

k=0 

k=1 

Q(i, j) 

m m-1  , 

= y A,(i,k)  cos  igi  j + y -F  5 

(5-12) 

k=0 

k=1 

These  equations  implicitly  define  the  cosine  and  sine  compo- 
nents U and  V of  P,  and  the  components  and  B-)  of  Q.  The 
reason  for  subscripting  A and  B will  become  apparent  later. 

Our  basic  procedure  will  be  to  derive  the  A's  and  B's  from  Q 
by  (5-12),  then  to  use  the  cyclic  reduction  process  to  obtain 
the  U's  and  V’s  from  the  A’s  and  B’s,  and  finally  to  synthesize 
P from  the  U’s  and  V’s  by  (5-11).  The  new  potential  distribu- 
tion P will  be  used  for  the  fields  in  the  next  trajectory  step, 
which  will  result  in  a new  space  charge  distribution  Q,  and  the 
process  is  repeated. 

Once  the  A^’s  and  B^’s  have  been  calculated,  they  contain  all 
the  information  about  the  space  charge  distribution  (in  a dif- 
ferent form) , and  matrix  Q can  be  vacated  and  used  for  storage 
of  the  U’s  and  V's  as  they  are  derived  from  the  A^’s  and  B^’s. 
Once  the  information  has  been  transferred  to  the  U’s  and  V’s, 
the  memory  space  for  the  A^'s  and  B'j’s  can  be  vacated,  and 
used  for  storage  of  the  potential  matrix  P.  Some  juggling  of 
indices  is  required,  but  there  is  a substantial  saving  of  mem- 
ory requirement.  The  combination  is 

U and  V share  Q’ s storage 

A^  and  B^  share  P ’ s storage . 
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The  orthogonality  of  the  Fourier  harmonics  allows  us  to  re- 
write (5-8)  as  2m  separate  equations;  there  are  m+1  cosine 
equations : 


(1+1/2i)  U(i+1,k)  cosTTkj/m  + (1-1/2i)  U(i-1,k)  cosTTkj/ra 
+ k^  U(i,k)|cosTtk(J-1)/m  + cosTtk(o+1  )/m| 


- 2(1+k^)  U(i,k)  cosTikj/m  = A^(i,k)  cosT^o/m  (5-13) 


for  k=0  to  m. 

There  are  m- 1 similar  sine  equations  relating  y and  , with 
sin  instead  of  cos,  running  k=1  to  m-1. 

But  by  a standard  trigonometric  identity 

cosTTk(j-1  )/m  + cosTTk(o+1)/m  = 2cosTTk3/m  cosT^/m  (5-14) 

with  a similar  sine  sm  formula.  Let  us  write 


F^(i,k)  = 1 - 1/2i  (5-15) 

G^(i,k)  = 1 + 1/2i  (5-16) 

S^(i,k)  = 2(1+f^)  - 2f^  cosTTk/m  ' (5-17) 

where  the  reason  for  the  dummy  k in  F>|  and  G-i , and  the  dummy 
i in  S'!  will  become  apparent  later. 

Substituting  (5-14)  in  (5-13),  simplifying  and  using  (5-15) 
through  (5-17)  we  have 

F^(i,k)  U(i-1,k)  - S^(i,k)  U(i,k)  + G^(i,k)  U(i+1,k)  = A^(i,k) 

(5-18] 

F^(i,k)  V(i-1,k)  - S^(i,k)  V(i,k)  + G.,(i,k)  V(i+1,k)  = B^(i,k) 

(5-19: 

where  (5-18)  runs  k=0  to  m,  (5-19)  runs  k=1  to  m-1,  and  both 
rian  i=1  to  Non-I. 


On  the  axis,  the  orthogonality  of  the  equi potentials  requires 


2U(1,k)  - S^(0,k)  U(0,k)  = A^(0,k),  k=0  to  m (5-20) 

2V(1,k)  - S^(0,k)  V(0,k)  = B^(0,k),  k=1  to  m-1  (5-21) 

5.3  Recursion  Step 

We  now  define  a sequence  of  functions  F^f  Gjjf  S^,  and 
by  the  following  recursion: 

Fji,k)  = F^_^(i,k).  F^_^(i-2'^-2,k)/S^_^(i-l'^"^,k)  (5-22) 

G^(i.k)  = G^_^(i,k).  G^_^(i+2'^-2,k)/S^_^(i+2’^-2,k)  (-25) 

S^(i,k)  = S^_^(i,k)  - F^_^(i,k).  G^_^(i-2'^"2,k)/S^_^(i-2'^-2,k) 

- G^_^(i,k).  F^_^(i+2^-2,k)/S^_^(1+2’^-2,k) 

(5-24) 

Aj^(i,k)  = + F^_^(i,k).  Aj^_^(i-2^-^,k)/S^_^(i-2^"2,k) 

+ G^_^(i,k).  Aj^_^(i+2^”^,k)/S^_^(i+2^"2,k) 

(5-25) 

\ 

+ Gjj_^(i,k).  Bjj_^(i+2"-2,k)/Sj^_^(i+2“-2,k) 

(5-26) 

The  recursion  runs  from  r=2  to  logo  (Nsr)  (i.e.  to  2,  3 or  4 
for  N3f^  = 4,  8 or  16),  for  i=2*^“"'  to  N3p^-2^^“"*  by  2^“%  and  for 
k=0  to  m.  (1  to  m-1  for  (5-26)).  The  n=1  values  have  already 
been  determined  by  (5-12)  and  (5-15)  through  (5-17);  the  reasons 
for  the  dummy  variables  in  tliese  should  now  be  clear. 
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5.4  Cyclic  Reduction 

We  now  write  down  (5-18)  for  consecutive  values  i-1,  i and  i+1 
(where  i is  even),  multiplying  the  first  by  F-^(i,k)/Si(i-1  ,k) 
and  the  third  by  (i,k)/S'](i+1  ,k) , and  adding.  We  apply 
(5-22)  through  (5-26)  v;ith  n=2.  Then  we  have,  after  simplifi- 
cation : 

F2(i,k)  U(i-2,k)  - 82(1, k)  U(i,k)  + G2(i,k)  U(i+2,k)  - A2(i,k) 

(5-27) 


Similarly  from  (5-19)  we  obtain 

F2(i,k)  V(i-2,k)  - S2(i,k)  V(i,k)  + G2(i,k)  V(i+2,k)  . B2(i,k) 

(5-28) 


Now  in  (5-27)  and  (5-28)  the  odd  numbered  rov/s  i-1  and  i+1 
have  been  eliminated,  and  the  form  of  the  equations  retained, 
with  n-2.  Thus  the  set  of  recursion  equations  (5-22)  through 
(5-26)  v;ere  properly  chosen,  and  can  be  applied  repeatedly  for 
successive  values  of  n.  The  general  forms  of  the  equations  are 


Fj^(i,k)  U(i-2^"\k)  - Sj^(i,k)  U(i,k)  + Gj^(i,k)  U(i+2’^"'' ,k) 

= Aj^(i,k)  (5-29) 

F^(i,k)  V(i-2^"\k)  - S^(i,k)  V(i,k)  + G^(i,k)  V(i+2’^-\k) 

= Bjj(i,k)  (5-30) 


running  i = 2^“"*  to  by  2^”"' 

k = 0 to  m for  (5-29),  1 to  m-1  for  (5-30)  . 

When  n = log2NsR  we  have  just  a single  pair  of  equations  for 
i = ^NgR,  which  we  will  write  ii/2  ior  brevity.  This  pair  of 
equations  relate  the  potentials  on  the  axis  and  at  the  wall  to 
the  row  ii/p  midway  along  the  tunnel  radius.  The  recursion 
stops  at  tnis  point,  but  we  still  do  not  know  the  potentials 
on  i'i/2  because  the  axial  potentials  are  not  known.  Taking 
the  wall  potential  to  be  zero,  the  final  pair  of  equations  is 
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The  resvJ-ting  simultaneous  solutions  are 


2U(2^,k)  - Tjj(k)  U(0,k)  = C^(k)  (5-39) 

2V(2^,k)  - T^(k)  V(0,k)  = Dj^(k)  (5-40) 

vhere  k = 0 to  m for  (5-33),  (5-34),  (5-36),  (5-37)  and  (5-39), 

= 1 to  m-1  for  (5-35),  (5-38)  and  (5-40). 

Again  the  recursion  stops  at  n = log^NsR,  and  we  have  (since 
U = V = 0 on  the  wall)  ^ 

U(0,k)  = -Cj^(k)/Tj^(k)  . (5-41) 

V(0,k)  = -Dj^(k)/T^(k)  (5-42) 


5.6  Backward  Recursion  and  Synthesis  of  P 


We  now  have  the  potential  components  on  the  axis;  these  can  be 
substituted  in  (5-31)  and  (5-32)  to  give  the  components  on  ii/2i 
the  mid  radius,  liien  the  values  on  the  axis  and  on  i^/?  deter- 
mine those  at  the  quarter  radius,  etc.  Specifically,  the  back- 
wards recursion  equations  are 


U(i,k)  = |Fj^(i,k)  U(i-2^^"\k)  - A^(i,k)  + Gj^(i,k)  U(i+2^"\k)|/ 

Sj^(i,k)  (5-43) 

V(i,k)  = (Fj.^(i,k)  V(i-2^"\k)  - Bj^(i,k)  + Gj^(i,k)  V(i+2^"\k)|/ 

Sj^(i,k)  (5-44) 


for  n = 1oS2^SR  tiy  -1 , 


n 

i. 

k 


2^"*’  to  Ngj^-2^"''  by  2^, 

0 to  m for  (5-43), 

1 to  m-1  for  (5-44)  . 
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On  completion  of  this  backwards  recursion,  we  have  all  the 
U and  V components.  We  now  insert  them  in  (5-11)  to  generate 
the  potential  matrix  P. 

While  these  recursions  appear  horrendous,  they  are  very  straight- 
forward for  computation.  Note  that  there  are  no  higher  functions 
to  be  evaluated,  because  the  sine  and  cosine  coefficients  in 
(5-11)  and  (5-12)  are  required  only  for  a fixed  set  of  submul- 
tiples of  2tt,  so  they  are  precalculated  and  stored  as  numerical 
coefficients  in  the  Fourier  analysis  and  synthesis  subroutines. 


5.7  Charge  Distribution  of  the  Beam 

The  foregoing  sections  5.1  to  5.6  have  shown  how  the  potential 
P is  derived  from  a given  charge  distribution  Q.  We  now  address 
the  question  of  exactly  what  distribution  best  represents  the 
beam.  The  test  of  what  is  ’best’  is  that  the  resulting  poten- 
tials, for  cases  to  which  an  analytic  solution  is  known,  should 
agree  with  the  analytic  values  as  closely  as  possible  for  as  wide 
a range  of  beam  diameters  as  possible.  Naturally  we  shall  find 
that  a finer  grid  — larger  values  of  Nsr  — will  give  better 
accuracy,  but  it  will  be  shown  that  accuracy  in  the  2%  region  can 
be  achieved  even  for  Nsf^=4  for  beams  v/ith  b/a  of  .5  or  greater, 
while  Nsr=8  gives  accuracy  better  than  1%  for  b/a  > .3,  which  is 
adequate  for  any  .foreseeable  TWT. 


5.8  The  Uniform  Beam 


We  consider  a uniform  beam  of  radius  b in  a tunnel  of  radius  a, 
at  voltage  Vq  and  microperveance  |jP;  the  potential  depression  on 
the  axis  is 


= -.0304  |.5  + tn(a/b)|,  (5-45) 

and  the  charge  density  (unifomi  out  to  r=b,  zero  for  r >b)  is 

p = 5.4x10'^°  V^(pP)/b  . (5-46) 

We  begin  with  a simple-minded  model,  in  which  each  node  of  the 
grid  lying  within  the  beam  (r  « b)  is  assigned  the  charge  Q 
given  by  (5-4)  and  (5-46)  combined,  and  each  node  outside  the 
beam  has  no  charge.  Obviously  this  model  will  give  errors  of 
one  sign  at  beam  diameters  such  that  a row  of  nodes  lies  just 
inside  the  beam,  and  of  the  opposite  sign  if  the  row  is  just 
outside.  Thus  we  expect  a sawtooth  curve  of  errors  as  a func- 
tion of  beam  radius.  Carrying  through  the  computation  for  a 
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case  with  Nc|^=l6,  we  obtain  a 'curve'  such  as  that  in  Figure  21. 

The  sawtooth  shape  is  not  quite  as  bad  as  it  looks:  the  errors 
are  within  10%  for  b/a  > .5»  and  there  would  be  a good  deal  of 
cancellation  of  positive  and  negative  errors  when  the  beam  diam- 
eter began  to  ch^ige,  as  it  will  in  the  real  tube. 

It  is  clear  that  the  curve  is  converging  to  a value  about  1% 
high  for  large  r/a.  This  is  a result  of  the  finite-difference 
treatment  of  the  problem.  The  offset  varies  as  it  is 

found  that  if  the  charges  assigned  to  the  nodes  are  all  reduced 
by  ttie  factor 

= 1 - 3.2/Ng^^  (5-47) 

then  the  offset  is  corrected  to  within  a fraction  of  1%  for 
Nsr  = 8,  16  or  32. 

At  the  left  of  Figure  21 , we  see  the  errors  becoming  quite  large 
for  small  b/a.  This  is  because  too  few  nodes  are  now  within  the 
beam  to  define  it  properly.  In  general,  we  find  that  the  errors 
for  large  b/a  depend  primarily  on  the  fineness  of  the  grid,  i.e. 
on  NsR,  while  the  errors  for  small  b/a  depend  mainly  on  the  ab- 
solute number  of  nodes  within  the  beam,  i.e.  on  (b/a)NsR. 

These  are  the  effect  of  the  finite  mesh  size;  when  we  bring  in 
the  discrete  super-electron  model  of  the  beam,  instead  of  the 
uniform  charge  density,  the  effects  are  more  complicated,  be- 
cause 'interference'  effects  arise  between  the  mesh  periodicity 
and  the  super-electron  periodicity.  The  super- electron  model  of 
the  beam  starts  with  a rectangular  array  of  super-electrons  in  a 
radial  plane.  The  number  of  columns  in  this  array  will  usually 
be  the  same  as  the  number  of  nodes  Nsa  (though  it  does  not  have 
to  be),  and  the  number  of  layers  Nl  will  normally  be  3 or  4, 
possibly  2,  6 or  even  8;  the  case  N3^=12,  Nsr=8,  Nl=4,  b/a=.7 
is  illustrated  in  Figure  22. 

If  we  simply  assign  the  charge  of  each  super-electron  to  its 
nearest  node,  we  shall  obtain  a sawtooth  error  curve  similar  to 
Figure  21.  To  avoid  this,  we  arrange  to  divide  the  charge  be- 
tween the  four  surrounding  nodes,  in  inverse  ratio  of  its  Ar 
and  Az  intercepts.  Referring  to  Figure  23,  if  the  super-elec- 
tron is  in  the  rectangle  defined  by  nodes  (i,j),  (i+1,j), 

(i,j+1)  and  (i+1,j+1)  (the  left  and  lower  sides  inclusive),  its 
charge  is  distributed  in  the  proportions 
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Figure  21 ; 


Errors  of  potential  on  axis  when  each  node 
inside  the  beam  is  assigned  the  full  charge 
given  by  (5-^),  and  each  node  outside  none. 


The  beam  is  assumed  to  be  uniform. 


16. 


Figure  22:  Superelectron  Distribution  for 

Uniform  Beam  (first  trial) 


Figure  23:  Distribution  of  Qiarge  of  a Super 
electron  to  Surrounding  Nodes 


(5 


Q(1-Az)(1-Ar) 

to 

(i,d) 

QAt(1- Az) 

to 

(i+1.a) 

QAz(1- Ar) 

to 

(i.:)+1) 

QAzAr 

to 

(i+1,j+1) 

where  Az  is  expressed  in  imits  of  the  axial  mesh  size  hgg^  and 
Ar  in  the  radial  mesh  size  hgr. 

Carrying  through  the  potential  calculation  for  this  case  we 
obtain  Figufe  24:  the  sawtooth  has  been  eliminated,  but  the 
offsets  at  the  wall  and  at  small  b/a  are  quite  large.  There 
is  nov/  only  a minor  dependence  on  Nsp^,  so  we  are  here  seeing 
mainly  the  effects  of  the  beam  model.  The  small  wrinkles  on 
the  curves  are  the  interference  effects  between  the  beam  and 
grid  models. 


Now  we  saw  in  the  continuous  charge  distribution  case  that  the 
offset  at  the  wall  could  be  corrected  by  slightly  adjusting  the 
amount  of  charge  assigned,  using  (5-46).  We  can  make  a similar 
type  of  correction  for  the  discrete  beam  model,  though  we  have 
not  been  able  to  find  an  analytic  expression  for  the  required 
reduction  factor,  which  we  will  call  q.  In  addition,  we  were 
clearly  incorrect  in  placing  the  outermost  layer  of  super-elec- 
trons in  Figure  22  ^ the  nominal  beam  radius.  If  they  are  to 
represent  the  real  electrons  in  their  neighborhood,  that  neigh- 
borhood should  sxirround  them;  thus  the  outermost  layer  should 
be  on  a line  at  some  radius  pb,  where  p < 1 (but  not  much  less), 
with  the  ethers  moved  in  proportionately. 


The  writer  has  no  analytic  method  of  finding  the  correction 
factors  p and  q,  but  computer  cut-and-try  is  effective.  It 
is  found  that  p affects  both  the  slope  and  the  absolute  level 
of  the  error  curves,  while  q mainly  affects  the  absolute  level. 
As  p is  reduced  from  unity,  the  slope  (over  the  interesting 
range  of  b/a)  decreases  and  eventually  changes  sign.  V/hen  a 
value  of  p is  fovind  which  makes  the  error  curve  as  flat  as 
possible,  q is  adjusted  to  level  it  around  zero.  Satisfactory 
(not  necessarily  optimum)  values  of  p and  q determined  in  this 
way  for  the  probable  combinations  of  Ncd  and  Nt,  are  given  in 
Table  5.1, 
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Figure  2k\  Errors  of  potential  on  the  axis  when  super- 
electron charges  are  distributed  to  the  4 
adjacent  nodes  instead  of  assignment  to  the 
nearest  one.  Beam  has  4 layers. 
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Figure  25:  Errors  of  potential  on  axis  when  correction 

factors  p and  q are  applied,  Ngn  = 2, 

2,  3 or  4 layers. 
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P E r c e n h 


Table  5.1  Correction  factors  for  super-electron  beam  model 
p is  correction  factor  for  radial  position,  q for 
charge 


P 

q 

— 

P 

q 

P 

q 

P 

q 

P 

16 

• 

.842 

.904 

.87"> 

.920 

.90 

.94 

.933 

.952 

8 

.770 

.889 

.855 

.916 

.870 

.934 

.90 

Co 

.917 

.955 

4 

.74 

.91 

.83 

.93 

.86 

.944 

.88 

.95 

.89 

.950 

2 

.65 

.90 

.72 

.90 

.75 

.905 

.78 

.95 

.80 

.920 

1 

2 

3 

4 

6 

8 J 

The  error  curves  obtained  using  these  factors  for  Nsr  =2,  4,  3 
and  16  are  shown  in  Figures  24,  25,  26  and  27.  The  last  is  good 
enough  for  any  conceivable  tube  design,  but  even  the  coarse  2 
mesh  case  (Figure  24)  v/ould  be  good  enough  for  preliminary  work 
with  b/a  > .5. 

Thus  Nsr  and  Nl  are  chosen  according  to  the  finenesj  of  model 
required,  remembering  that  both  directly  affect  the  computation 
time ; then  p and  q are  obtained  from  Table  5 . 1 » and  the  starting 
positions  for  the  super- electrons  are 


rg  = pbi/Nj^  i = 1 to  Nj^  (5-49) 

and  the  assigned  charges  are  (from  (5-4)) 

Qg  = qA^p  (b/a)  (Ngj^  Ksa/(Ni,€j,))  (5-50) 

= qfpbXg/(Nj_e^) 

where  p is  given  by  (5-46).  Note  that  the  correction  (5-4?) 
for  the  uniform  (fluid  model)' beam  is  not  applied  here,  as  it 
has  been  superseded  by  q. 


5.8.1  Accuracy  of  the  Uniform  Beam  Model 


From  Figures  24  through  2?  we  see  that  the  accuracy  for  the 
uniform  beam  can  be  expressed  by  the  following  tabulation: 


For  = 2 


errors  are  within  jf2?S  for  b/a  > .5 


4 


jf  1.496  for  b/a  > .45 


8 +196  for  b/a  > .33 

16  ±*596  for  b/a  > .22 


Clearly  the  accuracy  can  be  extended  to  smaller  beams  by 
using  larger  values  of  Nsr;  but  the  computation  time  will 
go  up  proportionately,  and  for  very  small  beams  a different 
approach  should  be  taken:  the  starting  point  should  be  a beam 
in  free  space,  v/ith  the  tunnel  introduced  only  as  a minor  per- 
turbation. fince  such  small  beams  are  not  of  interest  for  TWT 
work,  we  shall  not  pursue  this. 


5 . 9 The  Chopped  Beam  Accuracy  Check 

The  last  section  showed  that  very  satisfactory  accuracy  of 
the  potential  depression  is  obtained  for  the  uniform  beam, 
for  reasonable  values  of  Ncr;  but  this  did  not  check  the 
Fourier  analysis  part  of  the  procedure,  since  only  the  d-c 
term  remained. 

The  next  test  is  to  chop  the  beam  into  uniform  cylinders  of 
charge;  analytic  expression  for  the  potential  in  this  case 
are  given  by  Hechtel  [I6j,  RoWe  f17]  and  others.  If  the  disc 
thickness  is  equal  to  the  axial  mesh  size  b,  then  all  the  m 
Fourier  harmonics  are  required  to  express  the  potential,  and 
are  therefore  checked  for  accuracy. 


fl6]  J.R.  Hechtel,  'The  Effect  of  . otential  Beam  E*':orgy  on 

the  Performance  of  Linear  Beam  Devices',  lEEF  > ^’an  sections 
on  Electron  Devices,  ED-17,  i^11,  November  197  , pp.  999- 

1009. 

[17]]  J.E.  Rowe,  'Nonlinear  Electron-Wave  Interaction  Phenomena', 
Academic  Press,  New  York,  1965. 
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The  potentials  due  to  a single  disc  of  4 rings  located  at  one 
node  were  calculated  for  every  node  on  the  axis  (including  the 
zero  node  where  the  disc  was  located)  out  to  the  12th.  Beyond 
this  the  values  repeat,  of  course.  The  following  table  com- 
pares the  potentials  calculated  by  the  method  of  Sections  5.1 
through  5.6  with-  the  analytic  values  calculated  from  Hechtel's 
equations.  The  specific  case  chosen  was  a 30  kV  0.8  yP  beam 
for  which  the  analytic  values  Iiad  already  been  obtained. 

Nsa  =24,  Nsr=8,  Nl=4,  and  p and  q taken  from  Table  5.1. 


Table  5.2:  Potential  depression  on  axis  due  to  a 
single  disc  at  node  0 


ive  Node 

Pot.  dep. 
analytic  (volts) 

Pot.  Dep. 
FACE  (volts) 

0 

-148.40 

-146.25 

1 

-100.65 

-97.19 

2 

-57.72 

-57.98 

3 

-32.61 

-33.57 

4 

-18.33 

-19.24 

5 

-10.29 

-10.98 

6 

-5.77 

-6.26 

7 

-3.24 

-3.57 

3 

-1.83 

-2.04 

9 

-1.05 

-1.19 

10 

-.626 

-.722 

11 

-.419 

CM 

• 

12 

- .358 

-.422 

The  individual  discrepancies  are  nowhere  more  than  0.6%  of 
the  total  depression;  the  relative  discrepancies  are  somewhat 
higher,  but  this  is  largely  self-cancelling  when  the  whole  set 
of  24  discs  is  considered;  the  discrepancies  at  the  further 
nodes  are  quite  unimportant,  because  the  absolute  values  are 
here  so  small;  in  point  of  fact,  the  'analytic'  values  are 
suspect  here,  because  they  involve  calculation  of  a large 
number  of  Bessel  functions  which  nearly  cancel  each  other, 
so  that  round-off  errors  become  magnified.  The  total  of  all 
the  depressions  (including  the  min’ored  values  for  nodes  13 
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to  23)  should  a^ee  with  the  analytic  depression  (5-45)  of 
a uniform  beam  (since  24  adjacent  discs  constitute  such  a 
beam) ; we  find  that  the  FACR  total  actually  agrees  better 
with  (5-45)  than  the  'analytic'  total  of  Table  5.2. 

Thus  we  can  conclude  that  the  FACR  contributions  of  individ- 
ual discs  to  the  total  potential  are  accurate  to  better  than 
of  the  total  potential  depression,  and  are  at  least  as 
accurate  as  the  'analytic'  values,  for  Ng^=24. 

It  is  impossible,  of  course,  to  demonstrate  the  accuracy  of 
the  FACR  method  for  every  possible  case;  but  it  is  believed 
that  the  foregoing  checks,  for  a reasonably  typical  case, 
verify  that  the  method  is  sound,  and  that  no  mistakes  of 
scaling  have  been  made. 
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6.0  THE  TRAJECTORY  EQUATIONS 


The  foregoing  sections  have  provided  us  with  the  potential 
matrices,  and  a fast  and  accurate  interpolation  routine  by 
which  the  axial  and  radial  fields  can  be  derived  from  them 
at  the  position  of  each  superelectron  at  each  time  step. 
Thus  we  can  now  consider  the  fields  as  known. 

The  vector  equation  for  the  acceleration  of  an  electron  in 
combined  electric  and  magnetic  fields  is  quite  simple 


S = (e/m)E  + SXB  (6-1) 


In  principle,  one  can  separate  this  intc  three  component 
equations,  and  integrate  each  by  a Range-Kutta  or  similar 
routine.  This  procedure  is  inefficient  because  it  makes 
no  use  of  the  fact  that  we  know  jthe  integral  of  the^  equa- 
tion for  the  cross  field  case  (E  perpendicular  to  B) ; this 
is  the  well-known  cycloidal  solution,  combined  with  motion 
parallel  to  B which  is  not  affected  by  the  value  of  B. 
Textbook  formulations  of  the  cycloidal  solution  are  in  gen- 
eral too  simple  for  use  here  — they  do  not  allow  for  arbi- 
trary initial  velocity  components.  General  formulations  for 
the  cross  field  case  have  been  given  by  Yu  C'lS]  > Vaughan  [19] 
and  others.  For  the  present  purpose,  the  equations  in  [[19] 
are  the  more  convenient  starting  point,  since  they  give  the 
position  and  velocity  components  at  the  end  of  a time  step 
in  tenns  of  the  same  quantities  at  the  start,  together  with 
the  local  field  values. 

We  can  apply  that  formulation  to  the  present  case  by  adopt- 
ing a new  coordinate  system  as  shown  in  Figure  29.  Since  B 
has  no  o component  in  an  axisymmetric  system,  the  resultant 
B lies  in  the  r-z  plane,  and  the  new  Pz'  axis  is  taken  in 
this  direction;  Py'  is  normal  to  the  r-z  plane,  and  there- 
fox-e  makes  an  angle  0 with  the  Oy  direction,  and  Px'  is  then 


[18]  S.P.  Yu,  G.P.  Kooyers  and  0.  Buneman,  'Time-dependent  com- 
puter analysis  of  Electron-wave  Interaction  in  Crossed 
Fields',  J . Appl . Phy s . , vol.  36,  Aug.  1965,  pp  2550-2559. 

[19]  J.R.M.  Vaughan,  'Beam  Buildup  in  the  Dematron  Amplifier', 
IEEE  Transactions  on  Electron  Devices,  ED-18  #6,  June 
1971,  pp.  365-373. 


Figure  29:  Relation  of  Axial  and  Radial  Magnetic  Field 

Components  and  B„,  defining  the  auxiliay 
x' , y' , z’  coordinate  system.  lies  along 
the  resultant  B,  at  angle  0 to  r,  y'  is  nor- 
mal to  the  r-z  plane  at  P;  x’  in  the  r-z 
plane,  completing  the  triad. 


in  the  r-z  plane,  at  an  angle  0 (=  arc  tan  to  the 

radial  direction.  Thus  the  field  B is  normal  to  the  plane 
x'Py',  and  we  can  resolve  all  velocity  components  and  fields 
into  components  along  these  axes  and  use  the  solution  of  f19l. 
In  doing  so,  we. need  to  note  that  x'  and  y'  as  defined  here 
correspond  to  y and  x respectively  as  defined  there.  This 
formulation  is  extremely  accurate  — essentially  to  machine 
accuracy  — over  any  region  and  time  in  which  the  fields  are 
constant.  It  does  not  depend  on  At  being  small.  But  it  re- 
quires 57  multiplications,  31  additions  and  1 cosine  evalua- 
tion per  cycle.  The  complete  formulation  is  shown  in  Figtire 
30. 


This  is  more  than  we  need  for  the  present  purpose,  where  the 
fields  are  not  uniform  over  appreciable  distances  or  times. 

We  have  to  take  a At  that  is  a small  fraction  of  an  rf  per- 
iod; since  the  focusing  fields  are  always  such  that  me  is  of 
the  same  order  of  magnitude  as  m fv/ithin  a factor  of  about  3), 
we  shall  always  have  mcAt«1  also.  Hence  the  cos  (me  At)  and 
sin(u'cAt)  in  [192  can  be  replaced  by  1-a'c^At^/2  and  a'eAt  with 
extremely  small  errors.  These  aiid  consequent  simplifications 
were  carried  through  by  Prof.  0.  Buneman  in  1969,  for  the 
RZTRAJ  program  for  night  vision  devices  [3] , and  resulted  in 
a formulation  witi  only  31  multiplications',  21  additions  and 
no  cosine  e vaD.ua ti on  per  cycle  — a marked  and  valuable  reduc- 
tion, when  the  cycle  v;ill  be  repeated  several  thousand  times 
for  each  cavity  the  beam  is  tracked  through.  Extensive  tests 
reported  in  [3j  showed  that  this  compact  formulation  is  still 
capable  of  accuracy  at  the  or  better  level,  for  reasonable 

choices  of  step  size.  The  step  size  (At)  in  the  present  pro- 
gram is  not  specified  by  the  user,  but  is  automatically  chosen 
to  conform  to  the  accuracy  requirements  — it  is  tied  to  the 
number  of  radial  meshes  in  the  space  charge  matrix,  which  the 
user  can  specify.  If  the  user  sets  a large  value  (8  or  16) 
for  this,  indicating  that  he  desires  a high  accuracy  run  (and 
is  willing  to  pay  for  the  computer  time  involved),  a small  At 
will  result.  For  a 'debugging’  type  run  with  NgR  = 2,  a larger 
At  will  be  used,  but  still  not  large  enough  to  introduce  seri- 
ous errors  in  Buiieman’s  formulation. 

This  formulation,  which  is  the  one  used  in  the  program,  is 
shown  in  Figure  31 . It  has  eliminated  the  auxiliary  x' , y' , 
z'  coordinate  system,  but  it  is  so  tautly  written  that  without 
the  foregoing  tlieory  to  lead  up  to  it,  it  is  almoct  impossible 
to  see  how  it  works.  The  factor  A = 2/p+B1**2+B2**2+B3**2) 
is  the  approximation  to  the  missing  cosine  function,  with  a 
factor  2 resulting  from  dividing  the  step  into  two  parts. 


60b  Wl=SOH(Xr2+YI"‘2),  S2=YI /H I ,C2=X  l/Hl_  !_  RL^  SIN_PHI._COS  PHI 
610  B=SQK(133"2+Br2) , Cl  =33/B,SI  =B1/B  ! B TOT,  COS  THETA,  SIN  THETA 

620  00=E0*t3,A0=C5*[3  ! OMEGA  C = EB/M,  (F.3/M)(DT) 

62b  IF  A0<0.0002  THEN  I600__!  MAGNEFIC  F I hLD_ EFFECTIVELY  ZERO 

! 630  K6=C0S(A0)  ,K7=S0i<(  1 -K6*K6)  ,K0=  1 -K6  .Kh=!C7/B.K9=K0/B 
! 640  iCl=K7/OO.K2=KO/UO,K3=K2/B,K4  = {Kl-TO)/3 

i 64b  IF  Bl<.0001*B_niE;;  1600_  ! GUT0_ ROUTINE  FOR  AXIAL  FIELD  

j 6b0  o3=Sl*S2,C3=Cl *C2,S4=S2*CI ,C4=C2*S1 
1000  ! TRAJECTORIES 

1010  U5=U1*C3+V!  *S4-W1*S1 ^1  Ul' 

' 1020  Vb=Vl*C?-Ul *S2  ! VI' 

1030  W5=UUC4+Vl*S3+A'l  *C!  • ,i\' 

i 1100  E6=E3*C 1 +E1 >S 1 , Ej=-E3*S 1 +El *C| ! EZ' ,ER' 

1200  Xb=Ub*;<l -Vb*K2+Ej*K3  • X2' 

1210  Yb=U5*:<2+Vb*Kl-E-j*:<4  ! Y2' 

. 1220  Zb=^6*T0+E6*Kb  ! Z2'_ • ^ 

■ 1300  U6=Ub*K6+V5*K7  + Eo*i<8  ! U2' 

1310  V6=-*Ub*K/+Vb*;<6-L:b*K9  ! V2' 

; 1320  'A6=Hb+E6*Cb  J HZ'  _ 

1400  ! BACK  IJ  O.IIGINaL  COORDINATES 

1410  X2=X1 +Xb^C3-Yb*S2+Zb*C4 

I 1420  Y2=Y1 +Xb*S4+Y3*C2+Zb*S3_  _ 

I 1430  Z2=Zl+Zb*Cl -Xb^Si 
1500  U2=U6*C3-V6^S2+'R6*C4 

I 1510  V2=U6*J4+V6*C2+Ao*S3 

I 1520  W2=W6*C1-U6*S1 
! 1550  GOTO  2000 

1600  ! ROUriNb  FOR  Zl-.RO  MAGNETIC  FIELD 

1610  E2=Ei*K5 

1620  X2=X1+U1  ><T0+£2*C2 

. 1630  Y2=Y1+V1*T0-^E2*S2 ; 

1640  Z2=Z1+W1  xT0+E3*K-j 
1650  E2=E1*C5 

1660  U2=U1+E2*C2  

’ 1670  V2=V1+E2*S2 
1680  1V2=W1+E3«C5 

1690  GOTO  2000  _ _ „ 

IbOO  ! ROUTINE  FUR  AXIAL  MAGNET I'c  FIELD 
1820  X2=X1  +U  UKl  -VI  *K2-^E  1 *(K3*C2H(4*S2  ) 

1830  Y?=Y1+U1  *K2+V1*K1 +i;l  *(K3-^S2-K4*C2) . 

1840  Z2=Z1 +A’1*T0+E3*Kj 

1850  U2=U1*K6+V1 *K7+Ei *{K8*C2+K9*52) 

1660  V2=Vl*K6-Ul*K7+Et*(KoxrS2-K9*C2)_  . _ 

Ib/O  1\2=iM+E3*C5 

> 


Figure  30:  Full  Trajectory  Algorithm 
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33000  ? FIRST  STEP  INCREMENTS 

r 33200  El  =XO*Eb,E2=ro*t:5,E3=F2*T6  

35400  ? FIRST  INCREMENT  TO  STEP  SIZES 

35600  X3=X3+E1 ,Y3=Y3+u2,Z3=Z3+E3 

30  ' 36000  IF  ti7=l  THEN  39o00?  SIMPLIFIED  ROUTLNE  FOR  UNIFORM  FIELD 

37600  ? SCALED  MAGNETIC  FIELD  COMPONENTS 

3/bOO  bl=X0*E5,B2=Y0*£5,B3=F2*T6 
» 36000  ? CROSS  PRODUCTS  NITH  MAGNET IC_F I ELD 

36200  U=X3-Y3*i33+Z3*32,V=Y3-Z3*dl+X3*B3,R=Z3-X3*D2+Y3*i3l 
36400  B=2/(  I +Br2  + B2''2+B3''2) 

,?i  38600  IF  N>I  THEN  3y20C.  

36800  B=B/2 

39000  ? INCREMENT  WITH  CROSS  PRODUCTS 

_3l  39200  X3=X3-B*(V*B3-W*B2)  ,Y3=Y3-F*OV*Bl-U*B3),Z3=Z3-B*(U*i32-^Y*BI  ) 

39400  GUIO  40800 

39600  ? SIMPLIFIED  ROUTINE,  FOR  UNIFORM  FIELD 

« 39800  IF  B0=0  THEN  40800  _ 

40000  U=X3-Y3*B3,V=Y3fX3*B3 
40200  IF  N>l  THEN  40600 

il  40400  U=U/2,V=V/2  

40600  X3=X3-V*B4,Y3=Y3+U*B4 
40800  IF  N=l  THEN  41600 

« 41000  ? SECOND  HALF  STEP , FLECTRIC  FIELDS 

41200  X3=X3i-EI  ,Y3=Y3+u2,Z3=Z3+E3 
41400  ? INCREMENT  THE  COORDINATES 

il  41600  X0=X0+X3»Y0=Y0+V3,Z0=Z0+Z3/ZI=ZI +Z3  

41o00  RO=SQR(XO"2+YO''2  ) 


Figure  31:  Shortened  (Buneman)  Trajectory  Algorithm 


Under  no  circumstances  should  this  routine  be  altered.  It 
could  be  replaced  completely  by  the  more  general  routine  of 
Figure  30,  but  it  is  not  believed  that  there  would  be  any 
significant  increase  in  accuracy  of  the  overall  program  to 
compensate  for  the  shai-p  increase  in  running  time  that  must 
result. 


6. 1 Accuracy  of  the  Tra.jectory  Equations 

The  accuracy  of  Prof.  Buneman's  trajectory  algorithm  was 
tested  tiioroughly  in  the  development  of  the  image  intensi- 
fier  program  RZTRAJ,  and  the  tests  were  discussed  in  detail 
in  the  report  on  that  contract 

The  algorithm  could  not  be  taken  over  completely  ’onchanged 
from  RZTRAJ:  in  that  program  the  mesh  was  square,  and  the 
same  mesh  was  used  for  both  electrostatic  and  magnetic  fields. 
Since  this  gave  a unique  mesh  size  a?  a convenient  unit  of 
distance,  the  problem  was  normalized  to  that  distance.  In 
the  present  case,  we  have  three  meshes,  all  of  them  in  gen- 
eral rectangular,  and  independently  dimensioned.  As  a result 
we  decided  not  to  normalize  the  problem  to  any  one  mesh  size, 
but  to  retain  the  coordinates  in  MKS  units.  This  does  not 
change  the  basic  structure  of  Buneman's  algorithm,  but  it  re- 
quires changes  in  all  the  coefficients  used  in  converting  the 
gradients  derived  from  the  interpolation  routine  INTRA.  For 
example,  when  a mesh  is  allowed  to  be  rectangular  instead  of 
square,  two  distinct  coefficients  are  needed  for  the  axial 
and  radial  forces  instead  of  one  common  one. 

To  demonstrate  that  the  revised  routine  is  correct,  it  is 
primarily  necessary  to  show  that  the  expressions  used  for  the 
new  coefficients  are  correct;  if  they  are,  then  the  detailed 
checks  in  Q3J  will  apply  to  this  case  also. 

We  first  check  the  magnetic  field  coefficients  alone  by  track- 
ing a single  electron  in  a uniform  field,  so  that  the  path 
should  be  a circular  helix;  Figure  32  shows  the  projection  on 
the  X-Y  plane  of  such  a trajectory;  the  smaller  circle  results 
from  projecting  the  electron  with  the  transverse  velocity  com- 
ponent it  would  require  for  Brillouin  flow  if  the  entire  beam 
were  present.  In  the  absence  of  any  space  charge  force,  the 
diameter  of  the  circle  is  exactly  half  that  for  Brillouin  flow 
so  that  the  trajectory  passes  thiough  the  axis.  The  larger 
circle  results  from  doubling  the  transverse  injection  velocity 
Both  circles  are  tracked  round  about  ^00°,  and  in  the  overlap 
region  they  coincide  essentially  within  the  thickness  of  the 
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line.  This  demonstrates  that  no  mistake  has  been  made  in 
rescaling  tlie  magnetic  xield.  The  test  was  made  both  with 
the  ' uniform,  field-  bypass  in  the  trajectory  routine,  and 
with  a nominally  'nonuniform'  field  (so  that  the  full  tra- 
jectory routine  is  used),  but  the  field  was  in  fact  made 
uniform  to  a high  degree  by  specifying  a Helmholtz  pair  for 
the  coils.  The  differences  betv/een  the  two  cases  were  only 
at  the  0.1%  level,  comparable  to  the  differences  between  the 
Helmholtz  field  and  the  exactly  uniform  field,  so  t’\at  both 
versions  of  the  trajectory  routine  were  checked. 

The  smaller  circle  in  Figure  32  also  shows  that  the  program 
retains  RZTRAJ's  capability  to  follow  a trajectory  through 
the  axis  without  blowing  up;  as  in  the  earlier  case,  an 
abort  will  occur  if  the  trajectory  not  only  passes  through 
the  axis,"  but  also  has  one  of  the  calculated  points  on  it 
exactly  at  R=0.  The  probability  of  this  occurring  in  any 
actual  calculation  is  so  small  that  it  v;ould  be  a waste  of 
computer  time  to  put  a test  for  R=0  at  every  step  of  every 
ring,  and  go  to  a bypass  routine  if  it  occurred.  If  the 
problem  ever  did  occur,  it  could  be  eliminated  simply  by 
rerunning  the  case  with  a minute  change  of  some  input  var- 
iable. 

Having  established  that  the  magnetic  field  is  correctly 
scaled,  we  can  check  the  space  charge  field  by  rijmiing  a 
Brillouin  flow  case,  in  which  the  space  charge  force  should 
Just  balance  the  magnetic  and  inertial  forces.  This  can  be 
done  without  modifying  the  program,  simply  by  specifying  a 
uniform  magnetic  field  at  the  Erillouin  value,  together  with 
a drive  level  so  low  that  rf  fields  never  become  significant. 
In  the  'Navtest'  case  we  used  942  gauss,  and  10”9  watt  drive 
(95  dB  below  the  saturation  drive  of  3.2  watts).  Incidental- 
ly, this  demonstrates  the  very  large  dynamic  range  of  the 
program:  it  still  gives  a satisfactory  small  signal  calcula- 
tion at  this  level.  At  the  output  end,  the  rf  voltage  is 
still  below  1 volt,  and  the  beam  modulation  is  only  .0008, 
so  the  conditions  for  Brillouin  flow  are  substantially  met. 

Figure  33  shows  the  r-z  plane  projection  of  the  trajectories 
for  this  case. 

Finally,  the  rf  scaling  can  be  checked  by  using  a very  stiff 
beam,  so  ttat  the  ring  trajectories  are  tightly  constrained 
to  stay  close  to  their  starting  radii  by  the  large  magnetic 
field.  In  this  case  we  revert  to  normal  saturation  drive, 
and  postulate  a uniform  magnetic  field  ten  times  the  Brillouin 
value.  The  corresponding  cathode  immersion  is  98.5%,  and 
TRANS  assigns  the  appropriate  small  transverse  velocity  com- 
ponents. Under  these  conditions,  the  exit  bunching  and 
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velocities  of  the  rings  should  be  similar  to  the  exit  bunch- 
ing and  velocities  of  the  iiscs  as  given  by  NAV118B.  They 
■will  not  be  identical,  because  the  rings  are  at  different 
radii  from  the  used  as  the  center  of  force  for  the 

discs;  but  the  test  is  sufficient  to  show  whether  any  gross 
errors  such  as  omission  of  terms  or  sign  changes  have  oc- 
curred. Figure  3^  compares  the  exit  data  of  the  disc  and 
ring  models,  ar.d  shov/s  .that  they  are  satisfactorily  similar 
with  this  ultra-stiff  beam.  (When  it  was  first  tried,  the 
test  was  very  obviously  not  passed,  and  the  cause  v;as  found 
to  be  a 90°  phase  error  in  transferring  the  rf  fields  from 
NAV118  to  RNGTRJ.  Thus  the  test  is  sensitive  enough  to  be 
useful.)  We  have,  of  course,  already  discussed  the  accuracy 
of  the  potentials  themselves  in  Sections  3 through  5»  so 
that  correctness  of  the  coefficients  transfers  this  accuracy 
to  the  complete  calculation.  These  tests  do  not  completely 
eliminate  the  possibility  that  a coefficient  might  be  off  by 
a few  percent,  but  it  seems  very  unlikely:  the  terms  making 
up  the  coefficients  are  mostly  very  large  or  very  small  num- 
bers, so  that  the  derivations  can  be  expected  to  be  either 
correct,  or  off  by  orders  of  magnitude;  the  tests  have  shown 
that  the  latter  is  certainly  not  the  case. 
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Figure  34;  Comparison  of  Ultra-Stiff  Ring  Model  with 
Disc  Model  Exit  Velocities 
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7.0  ASSEMBLY  OF  THE  PROGRAM 


Subroutines  corresponding  to  the  foregoing  sections  were 
assembled  into  a complete  program  as  shown  by  the  block 
diagiam  in  Figure  35.  It  would  not  be  correct  to  say  that 
the  subroutines  v/ere  based  on  the  analyses,  because  in  some 
cases  the  subroutines  were  v/ritten  first  and  then  translated 
into  algebra  for  this  report.  It  has  to  be  recognized  that 
the  formalism  of  Fortran  or  Super  Basic  is  just  as  valid  for 
the  solution  of  problems  i\\  physics  and  engineering  as  the 
formalism  of  algebra,  and  may  well  be  the  best  one  to  use 
when  the  end  product  is  to  be  a program,  not  a textbook. 

In  such  cases,  translating  the  program  into  algebra  is 
really  only  conformity  to  what  will  soon  have  become  an 
outmoded  tradition.  The  reader  will  have  noticed  in  par- 
ticular that  Section  6 of  this  report  has  almost  no  equa- 
tions, only  copies  of  the  full  and  condensed  Super  Basic 
trajectory  algorithms:  no  algebra  translation  of  these  al- 
gorithms exists. 

The  entire  program  was  first  written  in  Su^er  Basic  to  allow 
each  section  to  be  tested  on  a time- sharing  system  with  a 
minimum  of  delay.  In  this  version,  each  second  level  sub- 
routine in  Figure  35  is  a self-contained  program,  and  the 
’Main  program’  is  the  user  at  the  keyboard,  calling  the  pro- 
grams in  the  correct  order.  The  third  level  subroutines  of 
Figure  35  are  then  GOSUBs  within  the  second  level  programs. 

jfipart  from  the  essentially  trivial  (but  agonizing)  problems 
of  tracking  down  bugs,  one  defect  of  principle  came  to  light, 
which  ought  to  have  been  foreseen:  the  procedure  of  moving 
beam  elements  forward  or  back  by  one  wavelength  v whenever  nec- 
essary to  keep  them  within  the  space  charge  wavelength  being 
tracked  has  been  used  in  all  prior  programs  of  this  type  known 
to  the  writer,  and  was  adopted  here  as  indicated  on  page  4. 

Biat  all  prior  programs  have  been  restricted  to  uniform  mag- 
netic fields.  If  this  procedure  is  adopted  with  a PPM  field, 
one  finds  that  the  PPM  period  is  not  equal  to  the  space  charge 
wavelength,  so  that  a ring  moved  up  or  back  one  wavelength 
suddenly  finds  itself  in  a different  magnetic  field,  for  which 
its  transverse  component  of  velocity  is  incorrect.  As  a re- 
sult, angular  momentum  is  not  conserved.  If  one  arbitrarily 
changes  the  tangential  velocity  to  match  the  local  field,  then 
energy  is  not  conserved. 

To  escape  from  this  dilemma,  we  adopted  the  valuable  concept 
of  the  Ddppelganger.  The  Doppelganger  (for  those  unfamiliar 
with  German  folklore)  is  a ghost  of  a living  individual  which 
haunts  that  ;'ndividual  but  occasionally  goes  elsewhere,  and 
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can  conveniently  be  blamed  for  any  mishap.  The  application 
of  this  to  the  present  problem  was  to  provide  each  ring  with 
a Dbppelganger  (vector  TW()  in  the  program).  The  rings  them- 
selves are  now  not  moved  up  or  back  when  they  leave  the 
tracked  wavelength:  only  the  I^ppelgangers  are  detached  to 

I stay  in  that  v/avelength  to  pick  up  the  appropriate  space 

charge  forces;  the  other  forces  are  still  evaluated  at  the 
true  positions  of  the  rings  and  thus  see  no  discontinuities. 

I Each  Ddppelganger  shares  all  its  host's  coordinates  and  veloc- 

ity components  except  for  the  Z coordinate,  so  that  only  one 
additional  coordinate  (stored  in  TW)  is  needed  for  each  ring. 


I 

i 


With  this  approach,  it  appears  that  both  angular  momentum  and 
energy  are  satisfactorily  conserved  in  both  PPM  and  uniform 
magrietic  fields.  This  is  visually  evident  in  the  X-Y  plots, 
which  show  substantially  radial  motion  at  the  exit  plane  for 
all  rings  which  reach  that  plane. 


I 


I 
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8.0  OUTPUT  OF  THE  PROGRAM 

Figure  36(a)  throu^  (o)  shows  the  complete  output  for  a 
’minimum'  case  corresponding  to  the  test  case  'Navtest' ; 

(a)  repeats  the  input  data  to  make  the  case  self -document- 
ing. (b)  is  a miscellaneous  collection  of  derived  quan- 
tities. (c)  is  the  matrix  RFMAT  and  rf  fields,  (d)  is 
the  matrix  MAGMAT  and.  magnetic  fields,  (e)  and  (f)  are 
the  disc  model  run through,  which  is  somewhat  beyond  satu- 
ration because  the  drive  level  (3.2  watts)  had  been  deter- 
mined for  the  12  disc  model,  and  the  6 disc  value  is  usually 
1 or  2 dB  different,  (g)  is  the  conversion  to  rings.  The 
varying  signs  in  the  YDOT  column  arise  because  some  rings 
are  in  one  half  of  the  magnet  period  and  some  in  the  other 
half,  where  the  field  is  reversed.  Each  starting  YDOT  is 
individually  assigned  to  match  the  local  field,  (h)  is 
the  ring  trajectory  printout,  giving  R and  Z for  each  ring 
every  20th  step,  (i)  is  the  tabulation  of  exit  coordinates 
and  velocity  components,  (j)  is  the  tabulation  of  intercep- 
tion on  the  ferrules,  both  by  ring  number  and  by  ferrule 
number,  follov/ed  by  some  statistics  on  ti;e  exiting  rings. 

The  trajectories  for  this  case  arc  plotted  in  Figures  37 
and  38.  An  over-large  plot  step  size  was  chosen,  so  that 
these  plots  are  somewhat  angular;  smoother  ones  can  be  dravm, 
but  require  more  'core. 

The  printout  for  a full  scale  96  ring  case  ;ls  identical  jn 
format  to  that  shown,  but  simply  more  massive  and  unsuitable 
for  reproduction. 


Interesting  sequences  of  cases  have  been  run:  with  a uniform 
field  ^ the  Brillouin  value,  end  no  rf  (10"®  watt  drive), 
the  oeam  remains  well-behaved  and  there  is  no  interception. 
When  3.2  watts  drive  is  applied,  but  nothing  else  changed, 
the  beam  'blows  up'  as  expected,  with  collection  of  20  out 
of  48  rings  spread  over  6 ferrules.  With  a PPM  field  peak- 
ing at  1.6  X Brillouin  (1530  gauss),  the  interception  is  re- 
duced to  8 rings  out  of  48,  and  3060  gauss  further  reduces 
Lhis  to  4 rings  on  two  ferrules. 


The  'hour  glass'  X-Y  plots  in  PPM  fields  are  different  from 
what  we  expected.  It  is  difficult  to  make  out  too  much  from 
the  completed  plot;  it  is  far  more  informative  to  watch  the 
plotter  drawing  them.  The  corresponding  plots  for  a Brill- 
ouin flow  case  are  arcs  of  circles,  as  expect'd. 
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VAUGHAN'S  RING  MODEL  TWT  PROGRAM  FORTRAN  VERSION  OF  MAH.  !976 


fTPiTINP  INPUT  FDR  RING  MODEL  PROGRAM 

, PART  1 t GENERAL*  • 

1 CASE  IDENTIFICATION 

UUUUUUUU 

NAVIESf 

1 LINEAR  UNIiS  (IN  OR  MM) 

1 PARI  2:  TUBE  PHYSICAL  DESCRIPTION* 
ITEM  1 TUNNEL  DIAMETER 

IN 

0.2000 

1 2 CAVITY  PERIOD 

3 gap  length 

, 4 TOTAL  ff  Or  CAViTIES 

0.3000 

0.1000 

' t NO.  OF  SEVER  CAVITIES 

1 SEVER  CAVITIES 

i 6 RING  CALC  START  aT  CAV  # 

2 

12  13 

24 

✓ 

1 7 FIELD  INTENSIFICATION  FACTOR 

2.6000 

t 

, PART  3!  COLD  TEST  DATA* 

b PHASE  VEL  (M/SEC/1E7) 
9 IMPEDANCE  (OHMS) 

10  LOSS  (DB/CAV) 

10.0000 

10.0000 

JD.iOOO- 

PART  4*  MAGNETIC  FIELD  DATA* 
i 11  UNIFORM  (1)  OR  PERIODIC  12) 

I 13  PERIODIC  LENGTH 

32  NUMBER  OF  COILS  (MAJC-'^fO) 
COIL  R 

0.6000 

4 

Z 

M 

1 0.1 4S0 

! 2 0.M80 

t 3 0,1 4B0 

-0.1600 
0.1  600 
0.4600 

-168.1000 
169.7500 
-! 69.7500 

4 0.1480 

PART  b*  RE  DATA* 

14  FREQUENCY  (GHZ) 

0.7600 

10.000 

: 68. 1000 

16  DRIVE  POWER  (WATTS) 

3.20 

1 

PART  6*  BEAM  DATA* 

16  BEAM  VOLTAGE  (kV) 
! 17  MICROPERVEAKCE 

18  NOMINAL  B/A 

38.000 

1.100 

0.700 

PART  V*  BEAM  MODtL* 

19  tf  OF  DISCS  PER  WAVELENGTH 

20  # OF  RINGS  PER  DISC 

6 

4 

i PART  8*  POTENTIAL  MESH  DIMENSIONS* 
21  SCMAT*  RADIAL.  AXIAL  

4_.  . 6 . _ 

22  REMAT*  RADIAL,  AXIAL 

23  MAGMAT*  RADIAL,  AXIAL 

4 12 

4 1 6 

Figure  36(a) : Test  Case  Input  Data  (first  page  of  output) 
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SPACE  CilAHbE  DENSITY  _7.l40l4D-03  CB/M**3 
MEAN  POTENTIAL  DEFNESoION  876.8  VOLTS 
POTENTIAL  DEPRESSION  ON  AXIS  1088.6  VOLTS 
MEAN  BEAM  VELOCIT  Y^  10 .790*1  E /_M 'SEC  _ ___ 


BEAM  CURRENT  AMP 

BEAM  WAVELENGTH  0.42b  IN 
BRILLOUIN  FIELD  941. 6_ GAUSS  _ _ 

TIME  STEP  FOR  RING  CALC  5.88525  PICD5EC 
TIME  STEP  FOR  DISC  CALC  16.66666  PICOSEC 


Figure  36(b);  Derived  Quantities 
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Figure  36(c):  Subroutine  LALALA  Output.  The  potential  and  field  tables  are  optional 

printout  items,  selected  by  IRFPRT  and  IFLDPR. 
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Figure  36(d):  Subroutine  MAGVA  Output.  The  complete  potential  matrix  has  one  more 

column  on  the  right,  identical  to  column  2,  for  v;hich  there  was  not 
room  on  the  printer  used  here.  The  potential  and  field  tables  are 
optional  printout  items,  belected  by  MPRT  and  MFPRT. 


NAVnO*  PRELIMINARY  DISC  MODEL  CALCULATION 


I 


MEAN  POT.  DEP.  -876.81  VOLTS,  BEAM  VEL.  1.0 

LAMBDA  E=  _ _ 10.79  MM;  BETA  E= 882_.3 

DISK  CHARGE  1 . 3bdObD-i 0 CB 

PLASMA  WVLGTH  60, o6  MM;  PLASMA  FREQ  1.1786 

___REDUCED_PLASMA  FREQ  6.7739 
POT.  DEP.  FROM  SUM  OF  DISKS  -879.31  VOLTS 
VOLTS  AT  NODES  0 THRU  3 DUE  TO  DISC  AT  MODE  0* 

-506.936  -152.51  1_  -28.62b  _rl0.097 

BETA(-1)=  628.3  GAMMA(-1)=  592.3  GAMMA*A= 

Ml  = 0.8972  M2  = 0.6629 

TOTAL  IMPEDANCL=  _648.  OHMS  _ _ _ 

PIERCE'S  C = 0.08123  SMALL  B = 0.901 

QC  = 0.35880  SMALL  D = 0.030 

PHASE  SHIFT  PER  CA7IT7  = 274.3  DEG  _(  1 .524PI) 

VOLTAGE  ATTENUATION  FACTOR  PER  CAV. =0.9886 
1V1=  202.54  rt3=  44.09  R5“  3.69 


1 .0789  7D  08  M/S 


1 .1 7860D  10  R/S 
6.77395D  09  R/S 


1 .5045 


J7  I 

C/ V 

. VOLTS 

POWER 

GAIN 

I-FUND 

..,HEL  „ 

_ CAV. LOSS 

EFFIC. 

F.N^BAL. 

I 

NC. 

R-F 

WAITS 

DB 

/I-DC 

PHASE 

PHASE 

WATTS 

PCNT. 

PCNT. 

I 

33 

I 

I 

64.40  . 

__3,200 

_O..CO. 

a.  0000. 

0.0 

_18Q...0„ 

0.073 

0.00- 

.1  00.00 

! 

2 

63,83 

3.  M4 

-0.08 

0.0020 

356.9 

184.4 

0.072 

-0.00 

99.97 

I 

3 

68.35 

3.605 

0.52 

0.0064 

348.5 

21  2.7 

0.082 

0.00 

99.96 

J3 j 

_4 

82.  19 

__5..213. 

2. 1 2_ 

0.0110. 

_339.8 

236.5_ 

__  .0.119 

0.00 

99.98 

5 

1 06.48 

8.749 

4.37 

0.0161 

334.4 

253.5 

0.  199 

0.00 

100.01 

6 

137.55 

14.599 

6.59 

0.0200 

332.4 

263.5 

0.332 

0.00 

100.04 

7 

177.59  _ 

24.336 

8.8L 

0.0262 

330.9, 

265.!- 

0.554, 

0.01 

100.04 

8 

225.59 

39.269 

10.89 

0.0320 

329.8 

266.  1 

0.894 

0.01 

1 00.06 

9 

275.37 

58.51  1 

12.62 

0.0342 

328.4 

263.9 

1.332 

0.02 

100.03 

39 

lO 

337.01 

87.638 

__14.38 

0,0435 

325.8 

259.5- 

1 - 995., 

— 0,03 

99,96 

I I 

431 .83 

M 3,89  0 

16.53 

0.0661 

321.5 

255.9 

3.275 

0.05 

1 00 . 00 

! 

I 2 

64.40 

3.200 

-0-00 

0.0883 

31  3.  1 

266. 1 

0.073 

-0.00 

99.78 

4? 

13 

64.40 

3.200 

_.-0.00 

0,1102 

JM.5 

270.5_ 

0.073 

--0 . 00 . 

99.76 

I 4 

265.71 

54.477 

12.3) 

0.1236 

321.9 

272.3 

1.240 

0.02 

99.75 

i 5 

448.45 

155. 1 76 

16.86 

0.1148 

324.9 

274.2 

3.532 

0.05 

99.74 

' 16 

598.21 

276.121 

19.36 

0.1007 

324.7 

269  ._4_ 

6.285 

0.09 

99.62 

I 7 

74 1 . 84 

424.638 

21  .23 

0.1056 

321.5 

25  7 . 6 

9.666 

0.14 

99.58 

18 

907.18 

635.  (X)  8 

22.98 

0.1286 

31  6,  1 

240.8 

14.455 

0.20 

99.70 

ii 

19 

1203.05 

1 116,765 

25.43 

0.2093 

309.3 

251  .2 

25. 421. 

0.36- 

.-99.60 

20 

• 1579.02 

1923.8j1 

2 7.79 

0.2520 

305.2 

257.8 

43.792 

0.62 

100.01 

2l 

2034.24 

3193.030 

29.99 

0.3082 

302.2 

260.2 

72.682 

1.03 

1 00.66 

22 

258  7.  10 

_5164.4j3 

_ 32.08  _D. 3802 

299.4 

260.^ 

J17.557 

J . 67 

101  .28 

23 

3221 .52 

8007.910 

33.98 

0.4:j92 

295.7 

256.3 

182.282 

2.59 

1 01 .39 

24 

4174.25 

1 3444.800 

36.23 

0.6862 

290.2 

253.0 

306.041 

4.  34 

101  .29 

54 

TRANSFER  Z COORDINATEo «. 

0.1/5209  0. 1 7604 1 0. 

TRANSFER  VELOCrilES  (*ll 

1 7671'! 
E-7)J 

0.178032  O". 

l'  8087  r 

0 .1 84  oo  f 

. 

ii 

lO. 

81^04  9. 

98402  _9 

.78672  _ 

J 0.23752  _J  1 

.20927 

-11.49444 

Figure  36(e):  Subroutine  NAV118,  disc  model  program,  first  part. 

The  last  ^ lines  are  the  disc  positions  and  velocities 
vhich  will  form  the  starting  data  for  the  ring  model 
calculation. 
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j CAV.  " VOLTS  POWEii  GAIN  I-FUND  ABS  HEL  CAV.LOSS  EFFIC.  EN.BAL. 
NO.  R-F  WAITS  DB  /I -DC  PHASE  PHASE  WATTS  PCNT.  PCNT. 


2b  5237.83  21168.920  38.21  0.8366  283.4  246.8  481.864  6.84  100.73 

26  6457.78  32178.300  40.02  1.0828273.9237.6  732.468  10.39  99.64 

27  _7380.38  42029.530_  41  , 18_L.  !221_262..L  222.2 956.709 ^1  3.57_  ,96.39 

28  8011.59  49526.040  41.900.9902254.4216.9  1127.350  15.99  94.34 

29  7740.13  46226. 69C  4 1 . 60  0 .6322  248. 3 170.9  1052.248  14.93  89.18 

' 30  _6922.96  _369ai.,150,  40.63-Q.5769_,245^_2  J L9..9 84.1 . 794  _J  1. 94  ,_.87.00 

’ T = 2.11667D-09  SEC,  NO  = 127  STEPS 

' TOTAL  CAVITY  LOSSES  = 5984.386  J'l __ 

' TOTAL  SEVER  POWER  = 147.1  W 

I ELECTRONIC  EFFIC.  = 13.92  PCT 

RESIDUAL  BEAM  K.E.=  _,_7,3..07_  PCI 

Z COORDINA'IES,  MM» 

223.977  222.350  222.077  221.465  227.539  226.398 

1 GRID  COORDI NATES s , 

3.550555  2.645588  2.494014  2.153943  5.531154  4.896686 

VELOC I TI ES  ( M/SEC/ 1 E 7 ) t 

, 8.658630  3.412133  8^.180522  J0..^908069 9-^9,96846._5..879843  

RELATIVE  ENERGIES: 

i 0.64399  0.60785  0.57484  1.02206  0.85844  0.67732 

' DISTANCE  = 217.59  MM _ 


Figure  Subroutine  NAV118  continued;  this  case  is  overdriven, 

since  the  pov/er  peaks  at  49.5  kW  at  cavity  28  and  has 
fallen  to  37  kW  at  cavity  30.  It  illustrates  the  fact 
that  the  program  does  have  overdrive  capability,  though 
the  accuracy,  indicated  by  the  Energy  Balance  (far 
right  column),  is  beginning  to  degrade  beyond  saturation. 


THANS*  SUBDIVISION  Or  DISCS  INTO  RINGS  

AND  COMPUTATION  OF  CHARGE- 1 NiJtiPFNUENT  MATRICES 
UNCORRECTEU  RIfiG  CHARGE  (EQU.b-bO)  81.28407 
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Figure  36(g):  Subroutine  TRANS  Output.  The  discs  of  Fig.  36(e) 

are  subdivided  into  4 rings  each.  The  Z and  ZDOT 
values  correspond  to  the  last  4 lines  of  Fig.  36(e), 
with  the  conversion  of  units  indicated  above  the 
table.  The  time  step  was  given  in  Fig.  36(b). 
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Figure  36(h):  Runnir.g  tabulation  of  R and  Z coordinates  only 


Figure  36(1 ):  Full  coordinates  and  velocity  components  of  rings  which  reach 

the  exit  plane  (mid- tunnel  following  the  last  cavity). 


2. TABULATION  OF  COLLECTEUJIINGS 
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Figure  36(j);  Summary  Tabulations 


Figure  37:  R-Z  Trajectory  Plot,  PPM  Focused  Case.  The  plot  is  unduly  angular 

because  a coarse  plot  interval  v/as  used. 


REFERENCES 


J.R.M,  Vaughan,  * CalcuHation  of  Coupled-Cavity  TWT  Per- 
formance', IEEE  Transactions  on  Electron  Devices,  ED-22 
#10,  October  1975,  pp.  830-'890. 

'Three-Dimensional  Evaluation  of  Energy  Extraction  in 
Output  Cavities  of  Klystron  Amplifiers' , H.  G.  Kosmahl 
and  L.  U.  Albers,  IEEE  Transactions  on  Electron  Devices, 
ED-20  #10,  Oct.  1973,  PP.  883-890. 

'Electron  Ray- Tracing  Program  for  Image  Intonsifiers' , 
Final  Report,  Contract  DAAK02-67-C-0182,  bv  J.R.M. 
Vaughan  and  0;  Buneman,  Sept.  1970. 

H.K.  Detweiler  and  J.E.  Rowe,  'Electron  Dynamics  and 
Energy  Conversion  in  0-Type  Linear  Beam  Devices'  in  . 
'Advances  in  Microwaves',  Vol.  6,  1971,  Academic  Press, 
p.  35.  The  pair  of  equations  (l4)  on  p.  39  do  not 
satisfy  Laplace. 

J.R.M.  Vaughan,  'Representation  of  Axisymmetric  Magnetic 
Fields  in  Computer  Progranis' , IEEE  Transactions  on  Elec- 
tron Devices,  ED-19  #2,  February  1972,  pp.  144-151. 

J.R.M.  Vaughan,  'Methods  of  Finding  the  Parameters  of 
Ideal  Current  Loops  for  Computer  Simulation  of  Magnetic 
Fields',  IEEE  Transactions  on  Electron  Devices,  ED-21 
#5,  May  1974,  pp.  310-312. 

J.  Jeans,  'The  Mathematical  Theory  of  Electricity  and 
Magnetism',  Cambridge  Uhiv.  Press,  5th  Ed.  1933,  p.  443. 


F.  Sterzer  and  W.W.  Siekanowicz,  'The  Design  of  Periodic 
Permanent  Magnets  for  Focusing  of  Electron  Beams',  RCA 
Review,  Vol.  18,  pp.  39-59,  Mar.  1957. 


H.G.  Kosmahl  and  G.M.  Branch,  'Generalized  Representation 
of  Electric  Fields  in  Interaction  Gaps  of  Klystrons  aiid 
Traveling  Wave  Tubes' , IEEE  Transactions  on  Electron 
Devices',  ED-20  #7,  July  1973,  PP.  621-629. 


C.C.  Wang,  'Electromagnetic  Field  Inside  a Cylinder  with 
a Gap',  Journal  of  Applied  Physics,  1^,  June  1945, 
pp.  351-366. 


N.W.  McLachlan,  'Bessel  Functions  for  Engineers',  Oxford 
Univ.  Press,  p.  158,  equations  22  and  24. 


-109- 


12.  Diophantus,  ' Arithmetica* , Univ.  of  Alexandria,  Egypt, 
ca.  320,  trans.  S.  Stevin,  pub.  Elsevier,  Leyden,  1634. 
(Newton  collection.  Bender  Library,  Stanford  University). 

13.  M.  Abramowitz  and  I.  Stegun,  'Handbook  of  Mathematical 
Functions',  N.B.S.  Washington,  D.C.  1964  or  Dover  Publi- 
cations, New  York,  1965,  page  811. 

14.  R.W.  Hockney,  'The  potential  Calciilation  and  some  Appli- 
cations', in  'Methods  in  Computationeil  Physics',  Ed.  B. 
Alder  et  al.  Academic  Press,  New  York  1970. 

15.  F.S.  Shaw,  'Relaxation  Metaods',  Dover  Publications, 

New  York,  1953* 

16.  iT.R.  Hechtel,  'The  Effec;t  of  Potential  Beam  Energy  on  the 
Performance  of  Linear  Beam  Devices',  IEEE  Transactions  on 
Electron  Devices,  ED-17,  #11,  November  1970,  pp.  999-1009. 

17.  J.E.  Rowe,  'Nonlinear  Electron-Wave  Interaction  Phenomena*, 
Academic  Press,  New  York,  1965. 

18.  S.P.  Yu,  G.P.  Kooyers  and  0.  Buneman,  ' Time- Dependent  Com- 
puter Analysis  of  Electron-wave  Interaction  in  Crossed 
Fields',  J.Appl.Phys. , vol.  36,  Aug.  1965,  pp  2550-2559. 

19.  J.R.M.  Vaughan,  'Beam  Buildup  in  the  Dematron  Amplifier', 
IEEE  Transactions  on  Electron  Devices,  ED-18  #6,  June 
1971,  pp.  365-373. 


-110- 


